#ifndef GITHASH_PP
#define GITHASH_PP "unknown"
#endif
#ifndef GITDATE_PP
#define GITDATE_PP "unknown"
#endif
#ifndef BUILDDATE_PP
#define BUILDDATE_PP "unknown"
#endif

!> \brief A collection of utility routines used throughout FDS.

MODULE COMP_FUNCTIONS

USE PRECISION_PARAMETERS
IMPLICIT NONE

CONTAINS


!> \brief Returns the wall clock time in seconds.

REAL(EB) FUNCTION CURRENT_TIME()
USE MPI
CURRENT_TIME = MPI_WTIME()
END FUNCTION CURRENT_TIME

!> \brief Return the current time and date in ISO_8601 format.
!> \param DATE A character string containing the date and time of day.

SUBROUTINE GET_DATE_ISO_8601(DATE)

INTEGER :: VALUES(8)
CHARACTER(LABEL_LENGTH), INTENT(OUT) :: DATE
INTEGER :: HOURS_TZ
INTEGER :: MINUTES_TZ

CALL DATE_AND_TIME(VALUES=VALUES)

HOURS_TZ = VALUES(4)/60
MINUTES_TZ = MOD(VALUES(4),60)

WRITE(DATE,'(I4.4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2,".",I3.3,SP,I3.2,":",SS,I2.2)') &
        VALUES(1), &
        VALUES(2), VALUES(3), VALUES(5), VALUES(6), VALUES(7), &
        VALUES(8), HOURS_TZ, ABS(MINUTES_TZ)

END SUBROUTINE GET_DATE_ISO_8601


!> \brief Return the current time and date
!> \param DATE A character string containing the date and time of day

SUBROUTINE GET_DATE(DATE)

INTEGER :: DATE_TIME(8)
CHARACTER(10) :: BIG_BEN(3),MONTH
CHARACTER(LABEL_LENGTH), INTENT(OUT) :: DATE

CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)

SELECT CASE(DATE_TIME(2))
   CASE(1)
      MONTH='January'
   CASE(2)
      MONTH='February'
   CASE(3)
      MONTH='March'
   CASE(4)
      MONTH='April'
   CASE(5)
      MONTH='May'
   CASE(6)
      MONTH='June'
   CASE(7)
      MONTH='July'
   CASE(8)
      MONTH='August'
   CASE(9)
      MONTH='September'
   CASE(10)
      MONTH='October'
   CASE(11)
      MONTH='November'
   CASE(12)
      MONTH='December'
END SELECT

WRITE(DATE,'(A,I3,A,I4,2X,I2.2,A,I2.2,A,I2.2)') TRIM(MONTH),DATE_TIME(3),', ',DATE_TIME(1), &
                                                DATE_TIME(5),':',DATE_TIME(6),':',DATE_TIME(7)
END SUBROUTINE GET_DATE


!> \brief Stop the code gracefully after writing a message.
!> \param MESSAGE Character string containing an explanation for shutting down.
!> \param PROCESS_0_ONLY If .TRUE., only MPI process 0 should write the message.

SUBROUTINE SHUTDOWN(MESSAGE,PROCESS_0_ONLY)

USE GLOBAL_CONSTANTS, ONLY: MYID,LU_ERR,CHID,SETUP_STOP,STOP_STATUS,N_MPI_PROCESSES
CHARACTER(*) MESSAGE
LOGICAL, INTENT(IN), OPTIONAL :: PROCESS_0_ONLY
CHARACTER(20) :: FMT

IF (PRESENT(PROCESS_0_ONLY)) THEN
   WRITE(FMT,'("(/A,A,I",I0,",A,A,A)")') CEILING(LOG10(REAL(N_MPI_PROCESSES+1,EB))) + 1
   IF (.NOT.PROCESS_0_ONLY) WRITE(LU_ERR,FMT) TRIM(MESSAGE),' (MPI Process:',MYID,', CHID: ',TRIM(CHID),')'
   IF (PROCESS_0_ONLY .AND. MYID==0) WRITE(LU_ERR,'(/A,A,A,A)') TRIM(MESSAGE),' (CHID: ',TRIM(CHID),')'
ELSE
   IF (MYID==0) WRITE(LU_ERR,'(/A,A,A,A)') TRIM(MESSAGE),' (CHID: ',TRIM(CHID),')'
ENDIF

STOP_STATUS = SETUP_STOP

END SUBROUTINE SHUTDOWN


!> \brief Return current system memory usage.
!>
!> Return the memory used by the given process at the time of the call. It only works on a Linux machine because
!> it searches for the system file called '/proc/PID/status' and reads the VmRSS value (kB).
!> The DEVC output QUANTITY 'RAM' outputs this value in units of MB.
!> Normally, this routine just returns 0 because it is non-standard Fortran and we have disconnected the GETPID call.
!> To invoke it, uncomment the GETPID statement, and the 'USE IFPORT' statement if compiling with Intel Fortran.

SUBROUTINE SYSTEM_MEM_USAGE(VALUE_RSS)

!USE IFPORT  ! Intel Fortran extension library. This is needed for GETPID.
INTEGER, INTENT(OUT) :: VALUE_RSS
CHARACTER(200):: FILENAME=' '
CHARACTER(80) :: LINE
CHARACTER(8)  :: PID_CHAR=' '
INTEGER :: PID
LOGICAL :: IFXST

VALUE_RSS=0  ! return zero if the memory file is not found

PID = 0
!PID = GETPID() ! GETPID is non-standard Fortran. Leave it commented out for the moment.

WRITE(PID_CHAR,'(I8)') PID
FILENAME = '/proc/'//TRIM(ADJUSTL(PID_CHAR))//'/status'

INQUIRE(FILE=FILENAME,EXIST=IFXST)
IF (.NOT.IFXST) then
   WRITE (*,*) 'WARNING: Memory system file does not exist'
  RETURN
ENDIF

OPEN(1000,FILE=FILENAME,ACTION='READ')
DO
   READ (1000,'(A)',END=120) LINE
   IF (LINE(1:6)=='VmRSS:') THEN
      READ (LINE(7:),*) VALUE_RSS
      EXIT
   ENDIF
ENDDO
120 CONTINUE
CLOSE(1000)

END SUBROUTINE SYSTEM_MEM_USAGE


!> \brief Read the name of the FDS input file, which is the first argument after the fds command itself.

SUBROUTINE GET_INPUT_FILE
USE GLOBAL_CONSTANTS, ONLY: FN_INPUT
IF (FN_INPUT=='null') CALL GET_COMMAND_ARGUMENT(1,FN_INPUT)
END SUBROUTINE GET_INPUT_FILE


!> \brief Assign a unique integer to be used as a file logical unit.

INTEGER FUNCTION GET_FILE_NUMBER()
USE GLOBAL_CONSTANTS, ONLY: MYID,FILE_COUNTER
FILE_COUNTER(MYID) = FILE_COUNTER(MYID) + 1
GET_FILE_NUMBER = FILE_COUNTER(MYID)
END FUNCTION GET_FILE_NUMBER


!> \brief Look through the FDS input file for the given namelist variable and then stop at that line.
!> \param NAME The four character namelist variable
!> \param LU Logical unit of the FDS input file
!> \param IOS Error code

SUBROUTINE CHECKREAD(NAME,LU,IOS)

USE GLOBAL_CONSTANTS, ONLY: INPUT_FILE_LINE_NUMBER,STOP_STATUS,SETUP_STOP,MYID,LU_ERR
INTEGER :: II
INTEGER, INTENT(OUT) :: IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(4), INTENT(IN) :: NAME
CHARACTER(80) TEXT
IOS = 1

READLOOP: DO
   READ(LU,'(A)',END=10) TEXT
   INPUT_FILE_LINE_NUMBER = INPUT_FILE_LINE_NUMBER + 1
   TLOOP: DO II=1,72
      IF (TEXT(II:II)/='&' .AND. TEXT(II:II)/=' ') EXIT TLOOP
      IF (TEXT(II:II)=='&') THEN
         IF (TEXT(II+1:II+4)==NAME) THEN
            IF (TEXT(II+5:II+5)==' ') THEN
               BACKSPACE(LU)
               IOS = 0
               EXIT READLOOP
            ELSE
               IF (MYID==0) WRITE(LU_ERR,'(/A,I0,A,A)') 'ERROR: Input line ',INPUT_FILE_LINE_NUMBER,&
                                                      ' is not formatted properly; NAMELIST: ',NAME
               STOP_STATUS = SETUP_STOP
               IOS = 1
               EXIT READLOOP
            ENDIF
         ELSE
            CYCLE READLOOP
         ENDIF
      ENDIF
   ENDDO TLOOP
ENDDO READLOOP

10 RETURN
END SUBROUTINE CHECKREAD


!> \brief Look for odd or illegal characters in the FDS input file.
!> \param LU Logical Unit of FDS input file.
!> \param IOS Error code, 0 means the bad text was found, 1 means it was not.
!> \param TEXT Illegal character found in FDS input file.

SUBROUTINE SCAN_INPUT_FILE(LU,IOS,TEXT)

INTEGER, INTENT(OUT) :: IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(80), INTENT(OUT) :: TEXT

IOS = 1

READLOOP: DO
   READ(LU,'(A)',END=10) TEXT
   IF (IACHAR(TEXT(1:1))==13 .AND. TEXT(2:2)=='&') THEN
      IOS = 0
      EXIT READLOOP
   ENDIF
ENDDO READLOOP

10 RETURN
END SUBROUTINE SCAN_INPUT_FILE


!> \brief Look for a certain TEXT string in the input file.
!> \param LU Logical Unit of the FDS input file
!> \param TEXT Character string to search for
!> \param FOUND T if the TEXT is found; F if it is not

SUBROUTINE SEARCH_INPUT_FILE(LU,TEXT,FOUND)

! Look for TEXT in the input file.

INTEGER, INTENT(IN) :: LU
LOGICAL, INTENT(OUT) :: FOUND
CHARACTER(*), INTENT(IN) :: TEXT
CHARACTER(180) :: LINE
INTEGER :: IND

FOUND = .FALSE.
REWIND(LU)

READLOOP: DO
   READ(LU,'(A)',END=10) LINE
   IND = INDEX(LINE,TEXT)
   IF (IND>0) THEN
      FOUND = .TRUE.
      EXIT READLOOP
   ENDIF
ENDDO READLOOP

10 RETURN
END SUBROUTINE SEARCH_INPUT_FILE


!> \brief Read through a comma-delimited output file and stop when the time, T, exceeds the last time read.
!> \param LU Logical Unit of the file to be appended
!> \param N_TEXT_LINES Number of header lines in the file
!> \param T The time at which to append the data

SUBROUTINE APPEND_FILE(LU,N_TEXT_LINES,T)

USE GLOBAL_CONSTANTS, ONLY: CLIP_RESTART_FILES
INTEGER, INTENT(IN) :: LU,N_TEXT_LINES
INTEGER :: ITER
REAL(EB) :: FILE_TIME
REAL(EB), INTENT(IN) :: T

READ_TEXT_LINES: DO ITER=1,N_TEXT_LINES
   READ(LU,*)
ENDDO READ_TEXT_LINES

READ_LOOP: DO
   READ(LU,*,END=10) FILE_TIME
   IF (FILE_TIME>T .AND. CLIP_RESTART_FILES) THEN
      BACKSPACE(LU)
      EXIT READ_LOOP
   ENDIF
ENDDO READ_LOOP

RETURN
10 BACKSPACE(LU)

END SUBROUTINE APPEND_FILE


!> \brief Reorder an input sextuple XB if needed.
!> \param XB User-specified real sextuplet.

SUBROUTINE CHECK_XB(XB)
REAL(EB) :: DUMMY,XB(6)
INTEGER  :: I
DO I=1,5,2
   IF (XB(I)>XB(I+1)) THEN
      DUMMY   = XB(I)
      XB(I)   = XB(I+1)
      XB(I+1) = DUMMY
   ENDIF
ENDDO
END SUBROUTINE CHECK_XB


!> \brief Change the units of the output quantity if it is an integrated quantity.
!> \param QUANTITY Name of output quantity.
!> \param UNITS Quantity units to be changed.
!> \param SPATIAL_STATISTIC Type of spatial integration.
!> \param TEMPORAL_STATISTIC Type of time integration.
!> \param LU_ERR Logical Unit of error output file.

SUBROUTINE CHANGE_UNITS(QUANTITY,UNITS,SPATIAL_STATISTIC,TEMPORAL_STATISTIC,LU_ERR)

USE GLOBAL_CONSTANTS, ONLY: MYID
CHARACTER(LABEL_LENGTH), INTENT(IN) :: QUANTITY,SPATIAL_STATISTIC,TEMPORAL_STATISTIC
INTEGER, INTENT(IN) :: LU_ERR
CHARACTER(LABEL_LENGTH), INTENT(INOUT) :: UNITS
CHARACTER(LABEL_LENGTH) :: NEW_UNITS
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
INTEGER :: I,UNIT_L

UNIT_L = LEN(TRIM(UNITS))
NEW_UNITS = UNITS
I = 1
SELECT CASE (SPATIAL_STATISTIC)
   CASE('VOLUME INTEGRAL')
      I = INDEX(UNITS,'/m3')
      IF (I/=0) WRITE(NEW_UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
      IF (TRIM(UNITS)=='1/s') THEN
         NEW_UNITS = 'm3/s'
         I = 1
      ENDIF
      IF (TRIM(UNITS)=='Pa') THEN
         NEW_UNITS = 'J'
         I = 1
      ENDIF
   CASE('MASS INTEGRAL')
      I = INDEX(UNITS,'/kg')
      IF (I/=0) WRITE(NEW_UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE('AREA INTEGRAL','SURFACE INTEGRAL')
      I = INDEX(UNITS,'/m2')
      IF (I/=0) WRITE(NEW_UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE ('VOLUME')
      NEW_UNITS = 'm3'
      I=1
   CASE ('AREA','SURFACE AREA')
      NEW_UNITS = 'm2'
      I=1
   CASE ('MINLOC X','MINLOC Y','MINLOC Z','MAXLOC X','MAXLOC Y','MAXLOC Z')
      NEW_UNITS = 'm'
      I=1
END SELECT

UNITS = NEW_UNITS

IF (I==0) THEN
   WRITE(MESSAGE,'(A,A,A,A)')  'WARNING: Problem with units compatibility of SPATIAL_STATISITIC ',TRIM(SPATIAL_STATISTIC), &
                               ' with the QUANTITY ',TRIM(QUANTITY)
   IF (MYID==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
ENDIF

I = 1

SELECT CASE (TEMPORAL_STATISTIC)
   CASE('TIME INTEGRAL')
      I = INDEX(UNITS,'/s')
      IF (I/=0) THEN
         WRITE(NEW_UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+2:UNIT_L)
      ELSE
         I = INDEX(UNITS,'kW')
         IF (I==1) WRITE(NEW_UNITS,'(A,A)') 'kJ',UNITS(3:UNIT_L)
         IF (I>1)  WRITE(NEW_UNITS,'(A,A,A)') UNITS(1:I-1),'kJ',UNITS(I+2:UNIT_L)
      ENDIF
      IF (QUANTITY=='FED') I=1
END SELECT

UNITS = NEW_UNITS

IF (I==0) THEN
   WRITE(MESSAGE,'(A,A,A,A)')  'WARNING: Problem with units compatibility of TEMPORAL_STATISTIC ',TRIM(TEMPORAL_STATISTIC), &
                               ' with the QUANTITY ',TRIM(QUANTITY)
   IF (MYID==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
ENDIF

END SUBROUTINE CHANGE_UNITS


END MODULE COMP_FUNCTIONS


!> \brief A collection of routines that manage memory

MODULE MEMORY_FUNCTIONS

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE COMP_FUNCTIONS, ONLY : SHUTDOWN
IMPLICIT NONE

CONTAINS

!> \brief Constructs an error message if there was an error during allocatiion of an array
!> \param CodeSect Character string containing the subroutine or function the allocation statement is in
!> \param VarName Character string containing the name of the array being allocated
!> \param IZERO Error value returned by the ALLOCATE statment where 0 is no error

SUBROUTINE ChkMemErr(CodeSect,VarName,IZERO)

! Memory checking routine

CHARACTER(*), INTENT(IN) :: CodeSect, VarName
INTEGER IZERO
CHARACTER(MESSAGE_LENGTH) MESSAGE

IF (IZERO==0) RETURN

WRITE(MESSAGE,'(4A)') 'ERROR: Memory allocation failed for ', TRIM(VarName),' in the routine ',TRIM(CodeSect)
CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)

END SUBROUTINE ChkMemErr


!> \brief Changes the allocation of an array with DIMENSION 1
!> \param P Original array
!> \param N1 Lower bound of new allocation
!> \param N2 Upper bound of new allocation

FUNCTION REALLOCATE(P,N1,N2)

! Resize the array P

REAL(EB), POINTER, DIMENSION(:) :: P, REALLOCATE
INTEGER, INTENT(IN) :: N1,N2
INTEGER :: NOLD, IERR
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
ALLOCATE(REALLOCATE(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE(N1:NOLD+N1-1) = P(N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P)

END FUNCTION REALLOCATE


!> \brief Changes the allocation of a string array
!> \param P Original array
!> \param CLEN Length of string
!> \param N1 Lower bound of new allocation
!> \param N2 Upper bound of new allocation

FUNCTION REALLOCATE_CHARACTER_ARRAY(P,CLEN,N1,N2)

! Resize the character array P

INTEGER, INTENT(IN) :: N1,N2,CLEN
CHARACTER(CLEN), POINTER, DIMENSION(:) :: P, REALLOCATE_CHARACTER_ARRAY
INTEGER :: NOLD, IERR
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
ALLOCATE(REALLOCATE_CHARACTER_ARRAY(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE_CHARACTER_ARRAY'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE_CHARACTER_ARRAY(N1:NOLD+N1-1) = P(N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P)

END FUNCTION REALLOCATE_CHARACTER_ARRAY


!> \brief Changes the allocation of an array with DIMENSION 2
!> \param P Original array
!> \param M1 Lower bound of first dimension of new allocation
!> \param M2 Upper bound of first dimension of new allocation
!> \param N1 Lower bound of second dimension of new allocation
!> \param N2 Upper bound of second dimension of new allocation

FUNCTION REALLOCATE2D(P,M1,M2,N1,N2)

REAL(EB), POINTER, DIMENSION(:,:) :: P, REALLOCATE2D
INTEGER, INTENT(IN) :: M1,M2,N1,N2
INTEGER :: MOLD,NOLD,IERR
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
ALLOCATE(REALLOCATE2D(M1:M2,N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE2D'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
MOLD = MIN(SIZE(P,DIM=1), M2-M1+1)
NOLD = MIN(SIZE(P,DIM=2), N2-N1+1)
REALLOCATE2D(M1:MOLD+M1-1,N1:NOLD+N1-1) = P(M1:MOLD+M1-1,N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P)

END FUNCTION REALLOCATE2D


!> \brief Changes the allocation of the array STRING by adding 100 new entries
!> \param NM Mesh number

SUBROUTINE RE_ALLOCATE_STRINGS(NM)

CHARACTER(MESH_STRING_LENGTH), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M=>NULL()

M=>MESHES(NM)
ALLOCATE(DUMMY(1:M%N_STRINGS))
DUMMY = M%STRING
DEALLOCATE(M%STRING)
ALLOCATE(M%STRING(M%N_STRINGS_MAX+100))
M%STRING(1:M%N_STRINGS) = DUMMY(1:M%N_STRINGS)
M%N_STRINGS_MAX = M%N_STRINGS_MAX+100
DEALLOCATE(DUMMY)

END SUBROUTINE RE_ALLOCATE_STRINGS


!> \brief Re-allocate the derived type array P_ZONE

SUBROUTINE REALLOCATE_P_ZONE(N1,N2)

INTEGER, INTENT(IN) :: N1,N2
TYPE (P_ZONE_TYPE), DIMENSION(:), ALLOCATABLE :: P_ZONE_DUMMY

ALLOCATE(P_ZONE_DUMMY(1:N2))
IF (ALLOCATED(P_ZONE)) THEN
   P_ZONE_DUMMY(1:N1) = P_ZONE(1:N1)
   DEALLOCATE(P_ZONE)
ENDIF
ALLOCATE(P_ZONE(1:N2))
P_ZONE(1:N2) = P_ZONE_DUMMY(1:N2)
DEALLOCATE(P_ZONE_DUMMY)

END SUBROUTINE REALLOCATE_P_ZONE


!> \brief Determines the size of the ONE_D_M_AND_E_XFER type structure for a surface type
!> \param SURF_INDEX Index of the surface type

SUBROUTINE COMPUTE_ONE_D_STORAGE_DIMENSIONS(SURF_INDEX)

USE GLOBAL_CONSTANTS, ONLY: N_TRACKED_SPECIES,N_LP_ARRAY_INDICES,N_SURFACE_DENSITY_SPECIES
INTEGER, INTENT(IN) :: SURF_INDEX
TYPE (SURFACE_TYPE), POINTER :: SF

SF => SURFACE(SURF_INDEX)

SF%ONE_D_REALS_ARRAY_SIZE( 1) = N_TRACKED_SPECIES            ! M_DOT_G_PP_ACTUAL
SF%ONE_D_REALS_ARRAY_SIZE( 2) = SF%N_MATL                    ! M_DOT_S_PP
SF%ONE_D_REALS_ARRAY_SIZE( 3) = N_TRACKED_SPECIES            ! M_DOT_G_PP_ADJUST
SF%ONE_D_REALS_ARRAY_SIZE( 4) = SF%N_CELLS_MAX+1             ! X
SF%ONE_D_REALS_ARRAY_SIZE( 5) = SF%N_CELLS_MAX+2             ! TMP
SF%ONE_D_REALS_ARRAY_SIZE( 6) = (SF%N_CELLS_MAX+2)*SF%N_MATL ! RHO
SF%ONE_D_REALS_ARRAY_SIZE( 7) = SF%N_LAYERS                  ! LAYER_THICKNESS
SF%ONE_D_REALS_ARRAY_SIZE( 8) = SF%NRA*SF%NSB                ! ILW
SF%ONE_D_REALS_ARRAY_SIZE( 9) = SF%NSB                       ! IL
SF%ONE_D_REALS_ARRAY_SIZE(10) = N_TRACKED_SPECIES            ! ZZ_G
SF%ONE_D_REALS_ARRAY_SIZE(11) = N_TRACKED_SPECIES            ! ZZ_F
SF%ONE_D_REALS_ARRAY_SIZE(12) = N_TRACKED_SPECIES            ! RHO_D_F
SF%ONE_D_REALS_ARRAY_SIZE(13) = N_TRACKED_SPECIES            ! RHO_D_DZDN_F
SF%ONE_D_REALS_ARRAY_SIZE(14) = N_LP_ARRAY_INDICES           ! A_LP_MPUA
SF%ONE_D_REALS_ARRAY_SIZE(15) = N_LP_ARRAY_INDICES           ! LP_CPUA
SF%ONE_D_REALS_ARRAY_SIZE(16) = N_LP_ARRAY_INDICES           ! LP_MPUA
SF%ONE_D_REALS_ARRAY_SIZE(17) = N_LP_ARRAY_INDICES           ! LP_TEMP
SF%ONE_D_REALS_ARRAY_SIZE(18) = SF%N_CELLS_MAX               ! RHO_C_S
SF%ONE_D_REALS_ARRAY_SIZE(19) = SF%N_CELLS_MAX+2             ! K_S
SF%ONE_D_REALS_ARRAY_SIZE(20) = N_SURFACE_DENSITY_SPECIES    ! AWM_AEROSOL
SF%ONE_D_REALS_ARRAY_SIZE(21) = (SF%N_CELLS_MAX+2)*SF%N_MATL ! RHO_DOT
SF%ONE_D_REALS_ARRAY_SIZE(22) = (SF%N_CELLS_MAX+2)*SF%N_SPEC ! RHO (gas)
! Update last array size index also to the sum below

SF%ONE_D_INTEGERS_ARRAY_SIZE(1) = SF%N_LAYERS           ! N_LAYER_CELLS

SF%N_ONE_D_STORAGE_REALS    = N_ONE_D_SCALAR_REALS    + SUM(SF%ONE_D_REALS_ARRAY_SIZE(1:22))
SF%N_ONE_D_STORAGE_INTEGERS = N_ONE_D_SCALAR_INTEGERS + SUM(SF%ONE_D_INTEGERS_ARRAY_SIZE(1:1))
SF%N_ONE_D_STORAGE_LOGICALS = N_ONE_D_SCALAR_LOGICALS

END SUBROUTINE COMPUTE_ONE_D_STORAGE_DIMENSIONS


!> \brief Determines the size of the WALL type structure for a surface type
!> \param SURF_INDEX Index of the surface type

SUBROUTINE COMPUTE_WALL_STORAGE_DIMENSIONS(SURF_INDEX)

INTEGER, INTENT(IN) :: SURF_INDEX
TYPE (SURFACE_TYPE), POINTER :: SF

SF => SURFACE(SURF_INDEX)

SF%N_WALL_STORAGE_REALS    = N_WALL_SCALAR_REALS    + SF%N_ONE_D_STORAGE_REALS
SF%N_WALL_STORAGE_INTEGERS = N_WALL_SCALAR_INTEGERS + SF%N_ONE_D_STORAGE_INTEGERS
SF%N_WALL_STORAGE_LOGICALS = N_WALL_SCALAR_LOGICALS + SF%N_ONE_D_STORAGE_LOGICALS

END SUBROUTINE COMPUTE_WALL_STORAGE_DIMENSIONS


!> \brief Determines the size of the CFACE type structure for a surface type
!> \param SURF_INDEX Index of the surface type

SUBROUTINE COMPUTE_CFACE_STORAGE_DIMENSIONS(SURF_INDEX)

INTEGER, INTENT(IN) :: SURF_INDEX
TYPE (SURFACE_TYPE), POINTER :: SF

SF => SURFACE(SURF_INDEX)

SF%N_CFACE_STORAGE_REALS    = N_CFACE_SCALAR_REALS    + SF%N_ONE_D_STORAGE_REALS
SF%N_CFACE_STORAGE_INTEGERS = N_CFACE_SCALAR_INTEGERS + SF%N_ONE_D_STORAGE_INTEGERS
SF%N_CFACE_STORAGE_LOGICALS = N_CFACE_SCALAR_LOGICALS + SF%N_ONE_D_STORAGE_LOGICALS

END SUBROUTINE COMPUTE_CFACE_STORAGE_DIMENSIONS


!> \brief Determines the size of the LAGRAGIAN_PARTICLE type structure for the input particle class
!> \param LPC_INDEX Index of the particle class

SUBROUTINE COMPUTE_PARTICLE_STORAGE_DIMENSIONS(LPC_INDEX)

INTEGER, INTENT(IN) :: LPC_INDEX
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE (SURFACE_TYPE), POINTER :: SF

LPC => LAGRANGIAN_PARTICLE_CLASS(LPC_INDEX)
SF  => SURFACE(LPC%SURF_INDEX)

LPC%N_STORAGE_REALS    = N_PARTICLE_SCALAR_REALS    + SF%N_ONE_D_STORAGE_REALS
LPC%N_STORAGE_INTEGERS = N_PARTICLE_SCALAR_INTEGERS + SF%N_ONE_D_STORAGE_INTEGERS
LPC%N_STORAGE_LOGICALS = N_PARTICLE_SCALAR_LOGICALS + SF%N_ONE_D_STORAGE_LOGICALS

END SUBROUTINE COMPUTE_PARTICLE_STORAGE_DIMENSIONS


!> \brief Allocates storage for data associated with a single WALL cell or LAGRANGIAN_PARTICLE
!> \param NM Index to the current mesh
!> \param SURF_INDEX Surface type. Used if WALL_INDEX or CFACE_INDEX is set
!> \param LPC_INDEX Lagrangian particle class type. Used if LP_INDEX is set.
!> \param WALL_INDEX Index to the wall cell being allocated. Optional (one of WALL_INDEX, CFACE_INDEX, or LP_INDEX)
!> \param CFACE_INDEX Index to the cut cell face being allocated. Optional (one of WALL_INDEX, CFACE_INDEX, or LP_INDEX)
!> \param LP_INDEX Index to the Lagrangian particle being allocated. Optional (one of WALL_INDEX, CFACE_INDEX, or LP_INDEX)
!> \param TAG Unique indentifier for a particle used since a particle can move between meshes and reneter a mesh.
!> \param NEW_TAG Flag indicating TAG is for a new particle.

SUBROUTINE ALLOCATE_STORAGE(NM,SURF_INDEX,LPC_INDEX,WALL_INDEX,CFACE_INDEX,LP_INDEX,TAG,NEW_TAG)

USE GLOBAL_CONSTANTS, ONLY: SIGMA,TMPA4,RPI,T_BEGIN,DIRICHLET,INITIALIZATION_PHASE,PERIODIC_TEST,TMPA
LOGICAL, INTENT(IN), OPTIONAL :: NEW_TAG
INTEGER, INTENT(IN):: NM,SURF_INDEX
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,LP_INDEX,LPC_INDEX,TAG
INTEGER :: N_NEW_STORAGE_SLOTS,STORAGE_INDEX,I1,I2,I,OLD_INDEX,STORAGE_INDEX_OLD,STORAGE_INDEX_OPEN,N_STORAGE_SLOTS_OLD,NAS
LOGICAL :: NEW_SLOT
TYPE (STORAGE_TYPE), POINTER :: OS
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE (LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE (LAGRANGIAN_PARTICLE_TYPE), ALLOCATABLE, DIMENSION(:) :: LP_DUMMY
TYPE (WALL_TYPE), ALLOCATABLE, DIMENSION(:) :: WALL_DUMMY
TYPE (CFACE_TYPE), ALLOCATABLE, DIMENSION(:) :: CFACE_DUMMY
TYPE (MESH_TYPE), POINTER :: M

M  => MESHES(NM)
SF => SURFACE(SURF_INDEX)

IF (PERIODIC_TEST==7) THEN
   N_NEW_STORAGE_SLOTS=2*(M%IBAR*M%JBAR + M%IBAR*M%KBAR + M%JBAR*M%KBAR)
ELSE
   N_NEW_STORAGE_SLOTS=1000
ENDIF

IF (PRESENT(WALL_INDEX)) THEN

   OS => MESHES(NM)%WALL_STORAGE(SURF_INDEX)

   ! Check if there are enough wall cells

   IF (WALL_INDEX>M%N_WALL_CELLS_DIM) THEN
      ALLOCATE(WALL_DUMMY(0:M%N_WALL_CELLS_DIM))
      WALL_DUMMY = M%WALL
      DEALLOCATE(M%WALL)
      ALLOCATE(M%WALL(0:M%N_WALL_CELLS_DIM+N_NEW_STORAGE_SLOTS))
      M%WALL(0:MESHES(NM)%N_WALL_CELLS_DIM) = WALL_DUMMY
      DEALLOCATE(WALL_DUMMY)
      M%N_WALL_CELLS_DIM = M%N_WALL_CELLS_DIM + N_NEW_STORAGE_SLOTS
   ENDIF

ELSEIF (PRESENT(CFACE_INDEX)) THEN

   OS => MESHES(NM)%CFACE_STORAGE(SURF_INDEX)

   ! Check if there are enough CFACEs

   IF (CFACE_INDEX>M%N_CFACE_CELLS_DIM) THEN
      ALLOCATE(CFACE_DUMMY(0:M%N_CFACE_CELLS_DIM))
      CFACE_DUMMY = M%CFACE
      DEALLOCATE(M%CFACE)
      ALLOCATE(M%CFACE(0:M%N_CFACE_CELLS_DIM+N_NEW_STORAGE_SLOTS))
      M%CFACE(0:MESHES(NM)%N_CFACE_CELLS_DIM) = CFACE_DUMMY
      DEALLOCATE(CFACE_DUMMY)
      M%N_CFACE_CELLS_DIM = M%N_CFACE_CELLS_DIM + N_NEW_STORAGE_SLOTS
   ENDIF

ELSEIF (PRESENT(LP_INDEX)) THEN

   LPC => LAGRANGIAN_PARTICLE_CLASS(LPC_INDEX)
   OS  => MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)

   ! Check if there are enough particles

   IF (LP_INDEX>MESHES(NM)%NLPDIM) THEN
      ALLOCATE(LP_DUMMY(MESHES(NM)%NLPDIM))
      LP_DUMMY = MESHES(NM)%LAGRANGIAN_PARTICLE
      DEALLOCATE(MESHES(NM)%LAGRANGIAN_PARTICLE)
      ALLOCATE(MESHES(NM)%LAGRANGIAN_PARTICLE(MESHES(NM)%NLPDIM+N_NEW_STORAGE_SLOTS))
      MESHES(NM)%LAGRANGIAN_PARTICLE(1:MESHES(NM)%NLPDIM) = LP_DUMMY
      DEALLOCATE(LP_DUMMY)
      MESHES(NM)%NLPDIM = MESHES(NM)%NLPDIM + N_NEW_STORAGE_SLOTS
   ENDIF

ELSE

   RETURN

ENDIF

! Determine if this is to be a new entry or old

STORAGE_INDEX_OLD  = 0
STORAGE_INDEX_OPEN = HUGE(STORAGE_INDEX_OPEN)

! If there is a specified particle TAG, find its storage slot

IF (PRESENT(TAG)) THEN
   IF (PRESENT(NEW_TAG) .AND. .NOT.NEW_TAG) THEN
      DO I=1,OS%N_STORAGE_SLOTS
         IF (OS%INTEGERS(1,I)==TAG) THEN
            STORAGE_INDEX_OLD = I
            EXIT
         ENDIF
      ENDDO
   ENDIF
ENDIF

! If there is an existing storage slot for the wall cell, use it

IF (PRESENT(WALL_INDEX) .AND. .NOT.INITIALIZATION_PHASE) THEN
   DO I=1,OS%N_STORAGE_SLOTS
      IF (OS%INTEGERS(1,I)==WALL_INDEX) THEN
         STORAGE_INDEX_OLD = I
         EXIT
      ENDIF
   ENDDO
ENDIF

! Look for the first available storage slot if it fails to find an existing one

IF (STORAGE_INDEX_OLD==0) THEN
   NAS = OS%NEXT_AVAILABLE_SLOT
   DO I=NAS,OS%N_STORAGE_SLOTS
      IF (OS%INTEGERS(1,I)==0) THEN
         STORAGE_INDEX_OPEN = MIN(I,STORAGE_INDEX_OPEN)
         OS%NEXT_AVAILABLE_SLOT = I+1
         EXIT
      ENDIF
   ENDDO
ENDIF

! Decide which storage index to use

IF (STORAGE_INDEX_OLD>0) THEN
   STORAGE_INDEX = STORAGE_INDEX_OLD
   NEW_SLOT = .FALSE.
ELSEIF (STORAGE_INDEX_OPEN<HUGE(STORAGE_INDEX_OPEN)) THEN
   STORAGE_INDEX = STORAGE_INDEX_OPEN
   NEW_SLOT = .TRUE.
ELSE
   STORAGE_INDEX = 0
   NEW_SLOT = .TRUE.
ENDIF

! Re-allocate storage arrays if necessary

IF (STORAGE_INDEX==0) THEN
   N_STORAGE_SLOTS_OLD = OS%N_STORAGE_SLOTS
   IF (PRESENT(LP_INDEX))    CALL REALLOCATE_STORAGE_ARRAYS(NM,1,LPC_INDEX,N_NEW_STORAGE_SLOTS)
   IF (PRESENT(WALL_INDEX))  CALL REALLOCATE_STORAGE_ARRAYS(NM,4,SURF_INDEX,N_NEW_STORAGE_SLOTS)
   IF (PRESENT(CFACE_INDEX)) CALL REALLOCATE_STORAGE_ARRAYS(NM,5,SURF_INDEX,N_NEW_STORAGE_SLOTS)
   DO I=1,N_STORAGE_SLOTS_OLD
      IF (PRESENT(WALL_INDEX))  OLD_INDEX = OS%INTEGERS(1,I)
      IF (PRESENT(CFACE_INDEX)) OLD_INDEX = OS%INTEGERS(1,I)
      IF (PRESENT(LP_INDEX))    OLD_INDEX = OS%INTEGERS(2,I)
      STORAGE_INDEX = I
      IF (PRESENT(WALL_INDEX)  .AND. OLD_INDEX>0) CALL WALL_POINTERS(OLD_INDEX,.FALSE.)
      IF (PRESENT(CFACE_INDEX) .AND. OLD_INDEX>0) CALL CFACE_POINTERS(OLD_INDEX,.FALSE.)
      IF (PRESENT(LP_INDEX)    .AND. OLD_INDEX>0) CALL PARTICLE_POINTERS(OLD_INDEX,.FALSE.)
   ENDDO
   STORAGE_INDEX = N_STORAGE_SLOTS_OLD + 1
ENDIF

! Assign values to the pointer arrays

IF (PRESENT(WALL_INDEX))  CALL WALL_POINTERS(WALL_INDEX,NEW_SLOT)
IF (PRESENT(CFACE_INDEX)) CALL CFACE_POINTERS(CFACE_INDEX,NEW_SLOT)
IF (PRESENT(LP_INDEX))    CALL PARTICLE_POINTERS(LP_INDEX,NEW_SLOT)

CONTAINS


!> \brief Assign pointers for variables associated with a LAGRANGIAN_PARTICLE
!> \param LP_INDEX_LOCAL Index of the LAGRANGIAN_PARTICLE in the array of particles for mesh NM
!> \param NEW Logical parameter indicating if the particle is new

SUBROUTINE PARTICLE_POINTERS(LP_INDEX_LOCAL,NEW)

INTEGER, INTENT(IN) :: LP_INDEX_LOCAL
LOGICAL, INTENT(IN) :: NEW

LP => MESHES(NM)%LAGRANGIAN_PARTICLE(LP_INDEX_LOCAL)

! Assign integer pointers and values

LP%TAG                 => OS%INTEGERS( 1,STORAGE_INDEX) ; IF (NEW) LP%TAG = TAG
LP%ARRAY_INDEX         => OS%INTEGERS( 2,STORAGE_INDEX) ;          LP%ARRAY_INDEX = LP_INDEX_LOCAL
LP%STORAGE_INDEX       => OS%INTEGERS( 3,STORAGE_INDEX) ; IF (NEW) LP%STORAGE_INDEX = STORAGE_INDEX
LP%CLASS_INDEX         => OS%INTEGERS( 4,STORAGE_INDEX) ; IF (NEW) LP%CLASS_INDEX = LPC_INDEX
LP%ORIENTATION_INDEX   => OS%INTEGERS( 5,STORAGE_INDEX) ; IF (NEW) LP%ORIENTATION_INDEX = 0
LP%WALL_INDEX          => OS%INTEGERS( 6,STORAGE_INDEX) ; IF (NEW) LP%WALL_INDEX = 0
LP%DUCT_INDEX          => OS%INTEGERS( 7,STORAGE_INDEX) ; IF (NEW) LP%DUCT_INDEX = 0
LP%DUCT_CELL_INDEX     => OS%INTEGERS( 8,STORAGE_INDEX) ; IF (NEW) LP%DUCT_CELL_INDEX = 0
LP%INIT_INDEX          => OS%INTEGERS( 9,STORAGE_INDEX) ; IF (NEW) LP%INIT_INDEX = 0
LP%CFACE_INDEX         => OS%INTEGERS(10,STORAGE_INDEX) ; IF (NEW) LP%CFACE_INDEX = 0

! Assign logical pointers and values

LP%SHOW              => OS%LOGICALS( 1,STORAGE_INDEX) ; IF (NEW) LP%SHOW  = .FALSE.
LP%SPLAT             => OS%LOGICALS( 2,STORAGE_INDEX) ; IF (NEW) LP%SPLAT = .FALSE.
LP%EMBER             => OS%LOGICALS( 3,STORAGE_INDEX) ; IF (NEW) LP%EMBER = .FALSE.
LP%PATH_PARTICLE     => OS%LOGICALS( 4,STORAGE_INDEX) ; IF (NEW) LP%PATH_PARTICLE = .FALSE.

! Assign real pointers and values

LP%X                 => OS%REALS( 1,STORAGE_INDEX) ; IF (NEW) LP%X       = 0._EB
LP%Y                 => OS%REALS( 2,STORAGE_INDEX) ; IF (NEW) LP%Y       = 0._EB
LP%Z                 => OS%REALS( 3,STORAGE_INDEX) ; IF (NEW) LP%Z       = 0._EB
LP%U                 => OS%REALS( 4,STORAGE_INDEX) ; IF (NEW) LP%U       = 0._EB
LP%V                 => OS%REALS( 5,STORAGE_INDEX) ; IF (NEW) LP%V       = 0._EB
LP%W                 => OS%REALS( 6,STORAGE_INDEX) ; IF (NEW) LP%W       = 0._EB
LP%PWT               => OS%REALS( 7,STORAGE_INDEX) ; IF (NEW) LP%PWT     = 1._EB
LP%ACCEL_X           => OS%REALS( 8,STORAGE_INDEX) ; IF (NEW) LP%ACCEL_X = 0._EB
LP%ACCEL_Y           => OS%REALS( 9,STORAGE_INDEX) ; IF (NEW) LP%ACCEL_Y = 0._EB
LP%ACCEL_Z           => OS%REALS(10,STORAGE_INDEX) ; IF (NEW) LP%ACCEL_Z = 0._EB
LP%RE                => OS%REALS(11,STORAGE_INDEX) ; IF (NEW) LP%RE      = 0._EB
LP%MASS              => OS%REALS(12,STORAGE_INDEX) ; IF (NEW) LP%MASS    = 0._EB
LP%T_INSERT          => OS%REALS(13,STORAGE_INDEX) ; IF (NEW) LP%T_INSERT= 0._EB
LP%DX                => OS%REALS(14,STORAGE_INDEX) ; IF (NEW) LP%DX      = 0._EB
LP%DY                => OS%REALS(15,STORAGE_INDEX) ; IF (NEW) LP%DY      = 0._EB
LP%DZ                => OS%REALS(16,STORAGE_INDEX) ; IF (NEW) LP%DZ      = 0._EB
LP%M_DOT             => OS%REALS(17,STORAGE_INDEX) ; IF (NEW) LP%M_DOT   = 0._EB
LP%HTC_LIMIT         => OS%REALS(18,STORAGE_INDEX) ; IF (NEW) LP%HTC_LIMIT=1.E6_EB

CALL ONE_D_POINTERS(LP_INDEX_LOCAL,NEW,N_PARTICLE_SCALAR_REALS,N_PARTICLE_SCALAR_INTEGERS,N_PARTICLE_SCALAR_LOGICALS)

END SUBROUTINE PARTICLE_POINTERS


!> \brief Assign pointers for variables associated with a WALL cell
!> \param WALL_INDEX_LOCAL Index of the WALL cell in the array of WALL cells for mesh NM
!> \param NEW Logical parameter indicating if the WALL cell is new

SUBROUTINE WALL_POINTERS(WALL_INDEX_LOCAL,NEW)

USE GLOBAL_CONSTANTS, ONLY: NULL_BOUNDARY,NEUMANN
INTEGER, INTENT(IN) :: WALL_INDEX_LOCAL
LOGICAL, INTENT(IN) :: NEW
TYPE (WALL_TYPE), POINTER :: WC

WC => MESHES(NM)%WALL(WALL_INDEX_LOCAL)

! Assign and initialize integers

WC%WALL_INDEX         => OS%INTEGERS( 1,STORAGE_INDEX) ; IF (NEW) WC%WALL_INDEX = WALL_INDEX
WC%SURF_INDEX         => OS%INTEGERS( 2,STORAGE_INDEX) ; IF (NEW) WC%SURF_INDEX = SURF_INDEX
WC%BACK_INDEX         => OS%INTEGERS( 3,STORAGE_INDEX) ; IF (NEW) WC%BACK_INDEX = 0
WC%BACK_MESH          => OS%INTEGERS( 4,STORAGE_INDEX) ; IF (NEW) WC%BACK_MESH  = NM
WC%BOUNDARY_TYPE      => OS%INTEGERS( 5,STORAGE_INDEX) ; IF (NEW) WC%BOUNDARY_TYPE = NULL_BOUNDARY
WC%SURF_INDEX_ORIG    => OS%INTEGERS( 6,STORAGE_INDEX) ; IF (NEW) WC%SURF_INDEX_ORIG = 0
WC%OBST_INDEX         => OS%INTEGERS( 7,STORAGE_INDEX) ; IF (NEW) WC%OBST_INDEX = 0
WC%PRESSURE_BC_INDEX  => OS%INTEGERS( 8,STORAGE_INDEX) ; IF (NEW) WC%PRESSURE_BC_INDEX = DIRICHLET
WC%LAPLACE_BC_INDEX   => OS%INTEGERS( 9,STORAGE_INDEX) ; IF (NEW) WC%LAPLACE_BC_INDEX = NEUMANN
WC%VENT_INDEX         => OS%INTEGERS(10,STORAGE_INDEX) ; IF (NEW) WC%VENT_INDEX = 0
WC%JD11_INDEX         => OS%INTEGERS(11,STORAGE_INDEX) ; IF (NEW) WC%JD11_INDEX = 0
WC%JD12_INDEX         => OS%INTEGERS(12,STORAGE_INDEX) ; IF (NEW) WC%JD12_INDEX = 0
WC%JD21_INDEX         => OS%INTEGERS(13,STORAGE_INDEX) ; IF (NEW) WC%JD21_INDEX = 0
WC%JD22_INDEX         => OS%INTEGERS(14,STORAGE_INDEX) ; IF (NEW) WC%JD22_INDEX = 0

! Assign and initialize reals

WC%DUNDT              => OS%REALS( 1,STORAGE_INDEX) ; IF (NEW) WC%DUNDT = 0._EB
WC%X                  => OS%REALS( 2,STORAGE_INDEX)
WC%Y                  => OS%REALS( 3,STORAGE_INDEX)
WC%Z                  => OS%REALS( 4,STORAGE_INDEX)
WC%VEL_ERR_NEW        => OS%REALS( 5,STORAGE_INDEX) ; IF (NEW) WC%VEL_ERR_NEW = 0._EB
WC%V_DEP              => OS%REALS( 6,STORAGE_INDEX) ; IF (NEW) WC%V_DEP = 0._EB
WC%Q_LEAK             => OS%REALS( 7,STORAGE_INDEX) ; IF (NEW) WC%Q_LEAK = 0._EB

CALL ONE_D_POINTERS(WALL_INDEX_LOCAL,NEW,N_WALL_SCALAR_REALS,N_WALL_SCALAR_INTEGERS,N_WALL_SCALAR_LOGICALS)

END SUBROUTINE WALL_POINTERS


!> \brief Assign pointers for variables associated with a CFACE
!> \param CFACE_INDEX_LOCAL Index of the CFACE in the array of CFACEs for mesh NM
!> \param NEW Logical parameter indicating if the CFACE is new

SUBROUTINE CFACE_POINTERS(CFACE_INDEX_LOCAL,NEW)

USE GLOBAL_CONSTANTS, ONLY: NULL_BOUNDARY
INTEGER, INTENT(IN) :: CFACE_INDEX_LOCAL
LOGICAL, INTENT(IN) :: NEW
TYPE (CFACE_TYPE), POINTER :: CFA

CFA => MESHES(NM)%CFACE(CFACE_INDEX_LOCAL)

! Assign and initialize integers

CFA%CFACE_INDEX        => OS%INTEGERS(1,STORAGE_INDEX) ; IF (NEW) CFA%CFACE_INDEX   = CFACE_INDEX
CFA%SURF_INDEX         => OS%INTEGERS(2,STORAGE_INDEX) ; IF (NEW) CFA%SURF_INDEX    = SURF_INDEX
CFA%VENT_INDEX         => OS%INTEGERS(3,STORAGE_INDEX) ; IF (NEW) CFA%VENT_INDEX    = 0
CFA%BACK_INDEX         => OS%INTEGERS(4,STORAGE_INDEX) ; IF (NEW) CFA%BACK_INDEX    = 0
CFA%BACK_MESH          => OS%INTEGERS(5,STORAGE_INDEX) ; IF (NEW) CFA%BACK_MESH     = NM
CFA%BOUNDARY_TYPE      => OS%INTEGERS(6,STORAGE_INDEX) ; IF (NEW) CFA%BOUNDARY_TYPE = NULL_BOUNDARY
CFA%CUT_FACE_IND1      => OS%INTEGERS(7,STORAGE_INDEX) ; IF (NEW) CFA%CUT_FACE_IND1 = -11
CFA%CUT_FACE_IND2      => OS%INTEGERS(8,STORAGE_INDEX) ; IF (NEW) CFA%CUT_FACE_IND2 = -11

! Assign and initialize reals

CFA%AREA               => OS%REALS( 1,STORAGE_INDEX)
CFA%X                  => OS%REALS( 2,STORAGE_INDEX)
CFA%Y                  => OS%REALS( 3,STORAGE_INDEX)
CFA%Z                  => OS%REALS( 4,STORAGE_INDEX)
CFA%NVEC(1:3)          => OS%REALS(5:7,STORAGE_INDEX)
CFA%VEL_ERR_NEW        => OS%REALS( 8,STORAGE_INDEX) ; IF (NEW) CFA%VEL_ERR_NEW = 0._EB
CFA%V_DEP              => OS%REALS( 9,STORAGE_INDEX) ; IF (NEW) CFA%V_DEP       = 0._EB
CFA%Q_LEAK             => OS%REALS(10,STORAGE_INDEX) ; IF (NEW) CFA%Q_LEAK      = 0._EB

CALL ONE_D_POINTERS(CFACE_INDEX_LOCAL,NEW,N_CFACE_SCALAR_REALS,N_CFACE_SCALAR_INTEGERS,N_CFACE_SCALAR_LOGICALS)

END SUBROUTINE CFACE_POINTERS


!> \brief Assign pointers for variables associated with the ONE_D derived type variable of a WALL, CFACE, or LAGRANGIAN_PARTICLE
!> \param ONE_D_INDEX_LOCAL Index of the WALL, CFACE, or LAGRANGIAN_PARTICLE
!> \param NEW Logical parameter indicating if the WALL, CFACE, or LAGRANGIAN_PARTICLE is new
!> \param RC Number of REALs already associated with the WALL, CFACE, or LAGRANGIAN_PARTICLE
!> \param IC Number of INTEGERs already associated with the WALL, CFACE, or LAGRANGIAN_PARTICLE
!> \param LC Number of LOGICALs already associated with the WALL, CFACE, or LAGRANGIAN_PARTICLE

SUBROUTINE ONE_D_POINTERS(ONE_D_INDEX_LOCAL,NEW,RC,IC,LC)

USE GLOBAL_CONSTANTS, ONLY: RADIATION,TMPM,RHOA,RSUM0,THERMALLY_THICK
INTEGER, INTENT(IN) :: ONE_D_INDEX_LOCAL,RC,IC,LC
LOGICAL, INTENT(IN) :: NEW
INTEGER :: N,NN
REAL(EB) :: X_POS
TYPE(ONE_D_M_AND_E_XFER_TYPE), POINTER :: ONE_D

IF (PRESENT(WALL_INDEX)) THEN
   ONE_D=>MESHES(NM)%WALL(ONE_D_INDEX_LOCAL)%ONE_D
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   ONE_D=>MESHES(NM)%CFACE(ONE_D_INDEX_LOCAL)%ONE_D
ELSEIF(PRESENT(LP_INDEX)) THEN
   ONE_D=>MESHES(NM)%LAGRANGIAN_PARTICLE(ONE_D_INDEX_LOCAL)%ONE_D
ENDIF

! Assign and initialize integers

ONE_D%ARRAY_INDEX     => OS%INTEGERS(IC+ 1,STORAGE_INDEX) ; IF (NEW) ONE_D%ARRAY_INDEX = ONE_D_INDEX_LOCAL
ONE_D%STORAGE_INDEX   => OS%INTEGERS(IC+ 2,STORAGE_INDEX) ; IF (NEW) ONE_D%STORAGE_INDEX = STORAGE_INDEX
ONE_D%II              => OS%INTEGERS(IC+ 3,STORAGE_INDEX) ; IF (NEW) ONE_D%II            = 0
ONE_D%JJ              => OS%INTEGERS(IC+ 4,STORAGE_INDEX) ; IF (NEW) ONE_D%JJ            = 0
ONE_D%KK              => OS%INTEGERS(IC+ 5,STORAGE_INDEX) ; IF (NEW) ONE_D%KK            = 0
ONE_D%IIG             => OS%INTEGERS(IC+ 6,STORAGE_INDEX) ; IF (NEW) ONE_D%IIG           = 0
ONE_D%JJG             => OS%INTEGERS(IC+ 7,STORAGE_INDEX) ; IF (NEW) ONE_D%JJG           = 0
ONE_D%KKG             => OS%INTEGERS(IC+ 8,STORAGE_INDEX) ; IF (NEW) ONE_D%KKG           = 0
ONE_D%IOR             => OS%INTEGERS(IC+ 9,STORAGE_INDEX) ; IF (NEW) ONE_D%IOR           = 0
ONE_D%PRESSURE_ZONE   => OS%INTEGERS(IC+10,STORAGE_INDEX) ; IF (NEW) ONE_D%PRESSURE_ZONE = 0
ONE_D%NODE_INDEX      => OS%INTEGERS(IC+11,STORAGE_INDEX) ; IF (NEW) ONE_D%NODE_INDEX    = 0
ONE_D%N_SUBSTEPS      => OS%INTEGERS(IC+12,STORAGE_INDEX) ; IF (NEW) ONE_D%N_SUBSTEPS    = 1

I1 = IC+1+N_ONE_D_SCALAR_INTEGERS ; I2 = I1 + SF%ONE_D_INTEGERS_ARRAY_SIZE(1) - 1
ONE_D%N_LAYER_CELLS(1:I2-I1+1) => OS%INTEGERS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS) = SF%N_LAYER_CELLS(1:SF%N_LAYERS)

! Assign and initialize logicals

ONE_D%BURNAWAY        => OS%LOGICALS(LC+1,STORAGE_INDEX) ; IF (NEW) ONE_D%BURNAWAY = .FALSE.

! Assign and initialize reals

ONE_D%AREA            => OS%REALS(RC+ 1,STORAGE_INDEX) ; IF (NEW) ONE_D%AREA            = 0._EB
ONE_D%HEAT_TRANS_COEF => OS%REALS(RC+ 2,STORAGE_INDEX) ; IF (NEW) ONE_D%HEAT_TRANS_COEF = 0._EB
ONE_D%Q_CON_F         => OS%REALS(RC+ 3,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_CON_F         = 0._EB
IF (RADIATION) THEN
   ONE_D%Q_RAD_IN     => OS%REALS(RC+ 4,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_RAD_IN        = SF%EMISSIVITY*SIGMA*TMPA4
   ONE_D%Q_RAD_OUT    => OS%REALS(RC+ 5,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_RAD_OUT       = SF%EMISSIVITY*SIGMA*TMPA4
ELSE
   ONE_D%Q_RAD_IN     => OS%REALS(RC+ 4,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_RAD_IN        = 0._EB
   ONE_D%Q_RAD_OUT    => OS%REALS(RC+ 5,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_RAD_OUT       = 0._EB
ENDIF
ONE_D%EMISSIVITY      => OS%REALS(RC+ 6,STORAGE_INDEX) ; IF (NEW) ONE_D%EMISSIVITY      = SF%EMISSIVITY
ONE_D%AREA_ADJUST     => OS%REALS(RC+ 7,STORAGE_INDEX) ; IF (NEW) ONE_D%AREA_ADJUST     = 1._EB
ONE_D%T_IGN           => OS%REALS(RC+ 8,STORAGE_INDEX) ; IF (NEW) ONE_D%T_IGN           = T_BEGIN
ONE_D%TMP_F           => OS%REALS(RC+ 9,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_F           = SF%TMP_FRONT
ONE_D%TMP_B           => OS%REALS(RC+10,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_B           = SF%TMP_INNER(SF%N_LAYERS)
ONE_D%U_NORMAL        => OS%REALS(RC+11,STORAGE_INDEX) ; IF (NEW) ONE_D%U_NORMAL        = 0._EB
ONE_D%U_NORMAL_S      => OS%REALS(RC+12,STORAGE_INDEX) ; IF (NEW) ONE_D%U_NORMAL_S      = 0._EB
ONE_D%U_NORMAL_0      => OS%REALS(RC+13,STORAGE_INDEX) ; IF (NEW) ONE_D%U_NORMAL_0      = 0._EB
ONE_D%RSUM_G          => OS%REALS(RC+14,STORAGE_INDEX) ; IF (NEW) ONE_D%RSUM_G          = RSUM0
ONE_D%TMP_G           => OS%REALS(RC+15,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_G           = TMPA
ONE_D%RHO_G           => OS%REALS(RC+16,STORAGE_INDEX) ; IF (NEW) ONE_D%RHO_G           = RHOA
ONE_D%U_TANG          => OS%REALS(RC+17,STORAGE_INDEX) ; IF (NEW) ONE_D%U_TANG          = 0._EB
ONE_D%RHO_F           => OS%REALS(RC+18,STORAGE_INDEX) ; IF (NEW) ONE_D%RHO_F           = RHOA
ONE_D%RDN             => OS%REALS(RC+19,STORAGE_INDEX) ; IF (NEW) ONE_D%RDN             = 1._EB
ONE_D%MU_G            => OS%REALS(RC+20,STORAGE_INDEX) ; IF (NEW) ONE_D%MU_G            = 0.1_EB
ONE_D%K_G             => OS%REALS(RC+21,STORAGE_INDEX) ; IF (NEW) ONE_D%K_G             = 0.1_EB
ONE_D%U_TAU           => OS%REALS(RC+22,STORAGE_INDEX) ; IF (NEW) ONE_D%U_TAU           = 0._EB
ONE_D%Y_PLUS          => OS%REALS(RC+23,STORAGE_INDEX) ; IF (NEW) ONE_D%Y_PLUS          = 1._EB
ONE_D%Z_STAR          => OS%REALS(RC+24,STORAGE_INDEX) ; IF (NEW) ONE_D%Z_STAR          = 1._EB
ONE_D%PHI_LS          => OS%REALS(RC+25,STORAGE_INDEX) ; IF (NEW) ONE_D%PHI_LS          = -1._EB
ONE_D%WORK1           => OS%REALS(RC+26,STORAGE_INDEX) ; IF (NEW) ONE_D%WORK1           = 0._EB
ONE_D%WORK2           => OS%REALS(RC+27,STORAGE_INDEX) ; IF (NEW) ONE_D%WORK2           = 0._EB
ONE_D%Q_DOT_G_PP      => OS%REALS(RC+28,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_DOT_G_PP      = 0._EB
ONE_D%Q_DOT_O2_PP     => OS%REALS(RC+29,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_DOT_O2_PP     = 0._EB
ONE_D%Q_CONDENSE      => OS%REALS(RC+30,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_CONDENSE      = 0._EB
ONE_D%TMP_F_OLD       => OS%REALS(RC+31,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_F_OLD       = SF%TMP_FRONT
ONE_D%K_SUPPRESSION   => OS%REALS(RC+32,STORAGE_INDEX) ; IF (NEW) ONE_D%K_SUPPRESSION   = 0._EB
ONE_D%BURN_DURATION   => OS%REALS(RC+33,STORAGE_INDEX) ; IF (NEW) ONE_D%BURN_DURATION   = SF%BURN_DURATION
ONE_D%T_SCALE         => OS%REALS(RC+34,STORAGE_INDEX) ; IF (NEW) ONE_D%T_SCALE         = 0._EB
ONE_D%Q_SCALE         => OS%REALS(RC+35,STORAGE_INDEX) ; IF (NEW) ONE_D%Q_SCALE         = 0._EB
ONE_D%L_OBUKHOV       => OS%REALS(RC+36,STORAGE_INDEX) ; IF (NEW) ONE_D%L_OBUKHOV       = 0._EB

I1 = RC+1+N_ONE_D_SCALAR_REALS ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(1) - 1
ONE_D%M_DOT_G_PP_ACTUAL(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%M_DOT_G_PP_ACTUAL(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(2) - 1
ONE_D%M_DOT_S_PP(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%M_DOT_S_PP(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(3) - 1
ONE_D%M_DOT_G_PP_ADJUST(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%M_DOT_G_PP_ADJUST(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(4) - 1
ONE_D%X(0:I2-I1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%X(0:SF%N_CELLS_INI) = SF%X_S(0:SF%N_CELLS_INI)

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(5) - 1
ONE_D%TMP(0:I2-I1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%TMP(0:SF%N_CELLS_INI+1) = SF%TMP_INNER(1)

IF (NEW) THEN
   IF (SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
      IF (SF%RAMP_T_I_INDEX > 0) THEN
         DO N=0,SF%N_CELLS_INI
            X_POS = MAX(0._EB,MIN(RAMPS(SF%RAMP_T_I_INDEX)%SPAN,SF%X_S(N)-RAMPS(SF%RAMP_T_I_INDEX)%T_MIN))
            ONE_D%TMP(N)=RAMPS(SF%RAMP_T_I_INDEX)%INTERPOLATED_DATA(NINT(X_POS*RAMPS(SF%RAMP_T_I_INDEX)%RDT)) + TMPM
         ENDDO
         ONE_D%TMP(SF%N_CELLS_INI+1) = ONE_D%TMP(SF%N_CELLS_INI)
         ONE_D%TMP_F = ONE_D%TMP(0)
         ONE_D%TMP_B = ONE_D%TMP(SF%N_CELLS_INI+1)
      ELSE
         ONE_D%TMP(0) = SF%TMP_INNER(1)
         ONE_D%TMP(SF%N_CELLS_INI+1) = SF%TMP_INNER(SF%N_LAYERS)
         NN=1
         DO N=1,SF%N_LAYERS
            ONE_D%TMP(NN:NN+SF%N_LAYER_CELLS(N)-1) = SF%TMP_INNER(N)
            NN = NN + SF%N_LAYER_CELLS(N)
         ENDDO
      ENDIF
   ELSE
      ONE_D%TMP(0:SF%N_CELLS_INI+1) = SF%TMP_INNER(1)
   ENDIF
ENDIF

IF (ALLOCATED(ONE_D%MATL_COMP)) DEALLOCATE(ONE_D%MATL_COMP)
ALLOCATE(ONE_D%MATL_COMP(1:SF%N_MATL))
DO NN=1,SF%N_MATL
   I1 = I2+1 ; I2 = I1 + (SF%N_CELLS_MAX+2) - 1
   ONE_D%MATL_COMP(NN)%RHO(0:SF%N_CELLS_MAX+1) => OS%REALS(I1:I2,STORAGE_INDEX)
   IF (NEW) ONE_D%MATL_COMP(NN)%RHO(0:SF%N_CELLS_INI+1) = SF%RHO_0(0:SF%N_CELLS_INI+1,NN)
ENDDO

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(7) - 1
ONE_D%LAYER_THICKNESS(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS) = SF%LAYER_THICKNESS(1:SF%N_LAYERS)

IF (ALLOCATED(ONE_D%BAND)) DEALLOCATE(ONE_D%BAND)
ALLOCATE(ONE_D%BAND(1:SF%NSB))
DO NN=1,SF%NSB
   I1 = I2+1 ; I2 = I1 + SF%NRA - 1
   ONE_D%BAND(NN)%ILW(1:SF%NRA) => OS%REALS(I1:I2,STORAGE_INDEX)
   IF (NEW) ONE_D%BAND(NN)%ILW(1:SF%NRA) = 0._EB
ENDDO

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(9) - 1
ONE_D%IL(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%IL(1:SF%NSB) = SIGMA*TMPA4/PI

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(10) - 1
ONE_D%ZZ_G(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%ZZ_G(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(11) - 1
ONE_D%ZZ_F(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%ZZ_F(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(12) - 1
ONE_D%RHO_D_F(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%RHO_D_F(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(13) - 1
ONE_D%RHO_D_DZDN_F(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%RHO_D_DZDN_F(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(14) - 1
ONE_D%A_LP_MPUA(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%A_LP_MPUA(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(15) - 1
ONE_D%LP_CPUA(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%LP_CPUA(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(16) - 1
ONE_D%LP_MPUA(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%LP_MPUA(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(17) - 1
ONE_D%LP_TEMP(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%LP_TEMP(1:I2-I1+1) = TMPA

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(18) - 1
ONE_D%RHO_C_S(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%RHO_C_S(1:I2-I1+1) = 1.E6_EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(19) - 1
ONE_D%K_S(0:SF%N_CELLS_INI+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%K_S(0:SF%N_CELLS_INI+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(20) - 1
ONE_D%AWM_AEROSOL(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%AWM_AEROSOL(1:I2-I1+1) = 0._EB

DO NN=1,SF%N_MATL
   I1 = I2+1 ; I2 = I1 + (SF%N_CELLS_MAX+2) - 1
   ONE_D%MATL_COMP(NN)%RHO_DOT(0:SF%N_CELLS_MAX+1) => OS%REALS(I1:I2,STORAGE_INDEX)
   IF (NEW) ONE_D%MATL_COMP(NN)%RHO_DOT(0:SF%N_CELLS_INI+1) = 0._EB
ENDDO

IF (ALLOCATED(ONE_D%SPEC_COMP)) DEALLOCATE(ONE_D%SPEC_COMP)
ALLOCATE(ONE_D%SPEC_COMP(1:SF%N_SPEC))
DO NN=1,SF%N_SPEC
   I1 = I2+1 ; I2 = I1 + (SF%N_CELLS_MAX+2) - 1
   ONE_D%SPEC_COMP(NN)%RHO(0:SF%N_CELLS_MAX+1) => OS%REALS(I1:I2,STORAGE_INDEX)
   IF (NEW) ONE_D%SPEC_COMP(NN)%RHO(0:SF%N_CELLS_INI+1) = 0._EB !Placeholder SF%RHO_0(0:SF%N_CELLS_INI+1,NN)
   !I1 = I2+1 ; I2 = I1 + (SF%N_CELLS_MAX+2) - 1
   !ONE_D%SPEC_COMP(NN)%RHO_DOT(0:SF%N_CELLS_MAX+1) => OS%REALS(I1:I2,STORAGE_INDEX)
   !IF (NEW) ONE_D%SPEC_COMP(NN)%RHO_DOT(0:SF%N_CELLS_INI+1) = 0._EB
ENDDO


END SUBROUTINE ONE_D_POINTERS

END SUBROUTINE ALLOCATE_STORAGE


!> \brief Increase the size of the array used to hold WALL, CFACE, LAGRANGIAN_PARTICLE data
!> \param NM Mesh number
!> \param MODE Indicator of the type of array
!> \param INDX Index of the LAGRANGIAN_PARTICLE_CLASS or the SURFACE type
!> \param N_NEW_STORAGE_SLOTS Number of new storage slots to allocate
!> \param NOM Optional parameter indicating the "number of the other mesh" for particles passing from one mesh to another

SUBROUTINE REALLOCATE_STORAGE_ARRAYS(NM,MODE,INDX,N_NEW_STORAGE_SLOTS,NOM)

INTEGER, INTENT(IN) :: MODE,NM,INDX,N_NEW_STORAGE_SLOTS
INTEGER, INTENT(IN), OPTIONAL :: NOM
INTEGER :: REALS_DIMENSION,INTEGERS_DIMENSION,LOGICALS_DIMENSION
TYPE (STORAGE_TYPE), POINTER :: OS
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DUMMY_REALS
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: DUMMY_INTEGERS
LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: DUMMY_LOGICALS

SELECT CASE(MODE)
   CASE(1:3)
      LPC => LAGRANGIAN_PARTICLE_CLASS(INDX)
      IF (MODE==1) OS  => MESHES(NM)%PARTICLE_STORAGE(INDX)
      IF (MODE==2) OS  => MESHES(NM)%OMESH(NOM)%ORPHAN_PARTICLE_STORAGE(INDX)
      IF (MODE==3) OS  => MESHES(NM)%OMESH(NOM)%ADOPT_PARTICLE_STORAGE(INDX)
      REALS_DIMENSION    = LPC%N_STORAGE_REALS
      INTEGERS_DIMENSION = LPC%N_STORAGE_INTEGERS
      LOGICALS_DIMENSION = LPC%N_STORAGE_LOGICALS
   CASE(4)
      SF => SURFACE(INDX)
      OS => MESHES(NM)%WALL_STORAGE(INDX)
      REALS_DIMENSION    = SF%N_WALL_STORAGE_REALS
      INTEGERS_DIMENSION = SF%N_WALL_STORAGE_INTEGERS
      LOGICALS_DIMENSION = SF%N_WALL_STORAGE_LOGICALS
   CASE(5)
      SF => SURFACE(INDX)
      OS => MESHES(NM)%CFACE_STORAGE(INDX)
      REALS_DIMENSION    = SF%N_CFACE_STORAGE_REALS
      INTEGERS_DIMENSION = SF%N_CFACE_STORAGE_INTEGERS
      LOGICALS_DIMENSION = SF%N_CFACE_STORAGE_LOGICALS
END SELECT

IF (OS%N_STORAGE_SLOTS>0) THEN
   ALLOCATE(DUMMY_REALS(REALS_DIMENSION,OS%N_STORAGE_SLOTS))
   ALLOCATE(DUMMY_INTEGERS(INTEGERS_DIMENSION,OS%N_STORAGE_SLOTS))
   ALLOCATE(DUMMY_LOGICALS(LOGICALS_DIMENSION,OS%N_STORAGE_SLOTS))
   DUMMY_REALS    = OS%REALS
   DUMMY_INTEGERS = OS%INTEGERS
   DUMMY_LOGICALS = OS%LOGICALS
   DEALLOCATE(OS%REALS)
   DEALLOCATE(OS%INTEGERS)
   DEALLOCATE(OS%LOGICALS)
ENDIF

ALLOCATE(OS%REALS(REALS_DIMENSION,OS%N_STORAGE_SLOTS+N_NEW_STORAGE_SLOTS))
ALLOCATE(OS%INTEGERS(INTEGERS_DIMENSION,OS%N_STORAGE_SLOTS+N_NEW_STORAGE_SLOTS))
ALLOCATE(OS%LOGICALS(LOGICALS_DIMENSION,OS%N_STORAGE_SLOTS+N_NEW_STORAGE_SLOTS))

OS%REALS    = 0._EB
OS%INTEGERS = 0
OS%LOGICALS = .FALSE.

IF (OS%N_STORAGE_SLOTS>0) THEN
   OS%REALS(1:REALS_DIMENSION,1:OS%N_STORAGE_SLOTS) = DUMMY_REALS
   OS%INTEGERS(1:INTEGERS_DIMENSION,1:OS%N_STORAGE_SLOTS) = DUMMY_INTEGERS
   OS%LOGICALS(1:LOGICALS_DIMENSION,1:OS%N_STORAGE_SLOTS) = DUMMY_LOGICALS
ENDIF

IF (OS%N_STORAGE_SLOTS>0) THEN
   DEALLOCATE(DUMMY_REALS)
   DEALLOCATE(DUMMY_INTEGERS)
   DEALLOCATE(DUMMY_LOGICALS)
ENDIF

OS%N_STORAGE_SLOTS = OS%N_STORAGE_SLOTS + N_NEW_STORAGE_SLOTS

END SUBROUTINE REALLOCATE_STORAGE_ARRAYS


!> \brief Find the index of the LAGRANGIAN_PARTICLE with the given TAG and LP_CLASS index. If not found, return 0.
!> \param NM Mesh number
!> \param LPC_INDEX Index of the LAGRANGIAN_PARTICLE_CLASS
!> \param LP_TAG A unique integer assigned to each LAGRANGIAN_PARTICLE, regardless of mesh
!> \param LP_INDEX The index of the LAGRANGIAN_PARTICLE of the input class and tag

SUBROUTINE GET_LAGRANGIAN_PARTICLE_INDEX(NM,LPC_INDEX,LP_TAG,LP_INDEX)

USE GLOBAL_CONSTANTS, ONLY: EVACUATION_ONLY
INTEGER, INTENT(IN) :: NM,LPC_INDEX,LP_TAG
INTEGER, INTENT(OUT) :: LP_INDEX
INTEGER :: I

LP_INDEX = 0
IF (EVACUATION_ONLY(NM)) RETURN
DO I=1,MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)%N_STORAGE_SLOTS
   IF (MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)%INTEGERS(1,I)==LP_TAG) THEN
      LP_INDEX = MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)%INTEGERS(2,I)
      EXIT
   ENDIF
ENDDO

END SUBROUTINE GET_LAGRANGIAN_PARTICLE_INDEX


END MODULE MEMORY_FUNCTIONS



!> \brief Functions and subroutines for manipulating geometry

MODULE GEOMETRY_FUNCTIONS

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
IMPLICIT NONE

CONTAINS


!> \brief Determine the mesh index and cell indices of an input point.
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param ZZ z-coordinate of the point (m)
!> \param NOM Number of the other mesh where the point is located
!> \param IIO I-index of the cell in the other mesh where the point is located
!> \param JJO J-index of the cell in the other mesh where the point is located
!> \param KKO K-index of the cell in the other mesh where the point is located

SUBROUTINE SEARCH_OTHER_MESHES(XX,YY,ZZ,NOM,IIO,JJO,KKO)

REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(OUT) :: NOM,IIO,JJO,KKO
TYPE (MESH_TYPE), POINTER :: M2=>NULL()

OTHER_MESH_LOOP: DO NOM=1,NMESHES
   IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP
   M2=>MESHES(NOM)
   IF (XX>=M2%XS .AND. XX<=M2%XF .AND.  YY>=M2%YS .AND. YY<=M2%YF .AND. ZZ>=M2%ZS .AND. ZZ<=M2%ZF) THEN
      IF (ALLOCATED(M2%CELLSI)) THEN
         XI  = MAX( 1._EB , MIN( REAL(M2%IBAR,EB)+ALMOST_ONE , M2%CELLSI(FLOOR((XX-M2%XS)*M2%RDXINT)) + 1._EB ) )
         YJ  = MAX( 1._EB , MIN( REAL(M2%JBAR,EB)+ALMOST_ONE , M2%CELLSJ(FLOOR((YY-M2%YS)*M2%RDYINT)) + 1._EB ) )
         ZK  = MAX( 1._EB , MIN( REAL(M2%KBAR,EB)+ALMOST_ONE , M2%CELLSK(FLOOR((ZZ-M2%ZS)*M2%RDZINT)) + 1._EB ) )
         IIO = FLOOR(XI)
         JJO = FLOOR(YJ)
         KKO = FLOOR(ZK)
      ELSE
         IIO = 0  ! The mesh if found, but no detailed information is available to the current process
         JJO = 0
         KKO = 0
      ENDIF
      RETURN
   ENDIF
ENDDO OTHER_MESH_LOOP

NOM = 0

END SUBROUTINE SEARCH_OTHER_MESHES


!> \brief Block or unblock a cell of a given mesh
!> \param NM Mesh number
!> \param I1 Lower bound of I cell indices of the region to be blocked or unblocked
!> \param I2 Upper bound of I cell indices of the region to be blocked or unblocked
!> \param J1 Lower bound of J cell indices of the region to be blocked or unblocked
!> \param J2 Upper bound of J cell indices of the region to be blocked or unblocked
!> \param K1 Lower bound of K cell indices of the region to be blocked or unblocked
!> \param K2 Upper bound of K cell indices of the region to be blocked or unblocked
!> \param IVAL Indicator if the cell is to blocked (1) or unblocked (0)
!> \param OBST_INDEX Index of the OBSTruction blocking the cell

SUBROUTINE BLOCK_CELL(NM,I1,I2,J1,J2,K1,K2,IVAL,OBST_INDEX)

INTEGER :: NM,I1,I2,J1,J2,K1,K2,IVAL,I,J,K,OBST_INDEX,IC
TYPE (MESH_TYPE), POINTER :: M=>NULL()

M => MESHES(NM)
DO K=K1,K2
   DO J=J1,J2
      DO I=I1,I2
         IC = M%CELL_INDEX(I,J,K)
         SELECT CASE(IVAL)
            CASE(0)
               M%SOLID(IC)        = .FALSE.
               M%OBST_INDEX_C(IC) = 0
            CASE(1)
               M%SOLID(IC)        = .TRUE.
               M%OBST_INDEX_C(IC) = OBST_INDEX
         END SELECT
         IF (OBST_INDEX==0) M%EXTERIOR(IC) = .TRUE.
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE BLOCK_CELL


!> \brief For each cell in an HT3D solid obstruction, determine its nearest surface WALL cell
!> \param NM Mesh number
!> \details Fill in the array MESHES(NM)\%WALL_INDEX_HT3D(IC,-3:3) such that: \n
!> WALL_INDEX_HT3D(IC,-3) is the wall index of the bottom of the solid obstruction \n
!> WALL_INDEX_HT3D(IC,-2) is the wall index of the back   of the solid obstruction \n
!> WALL_INDEX_HT3D(IC,-1) is the wall index of the left   of the solid obstruction \n
!> WALL_INDEX_HT3D(IC, 1) is the wall index of the right  of the solid obstruction \n
!> WALL_INDEX_HT3D(IC, 2) is the wall index of the front  of the solid obstruction \n
!> WALL_INDEX_HT3D(IC, 3) is the wall index of the top    of the solid obstruction \n
!> WALL_INDEX_HT3D(IC, 0) is the wall index of the nearest surface of the solid obstruction

SUBROUTINE ASSIGN_HT3D_WALL_INDICES(NM)

INTEGER, INTENT(IN) :: NM
INTEGER :: I,J,K,N,IC,II,JJ,KK,ICN,CELL_COUNT(-3:3)
TYPE(MESH_TYPE), POINTER :: M
TYPE(OBSTRUCTION_TYPE), POINTER :: OB

M => MESHES(NM)

OBST_LOOP: DO N=1,M%N_OBST
   OB => M%OBSTRUCTION(N)
   HT3D_IF: IF (OB%HT3D) THEN
      K_LOOP: DO K=OB%K1+1,OB%K2
         J_LOOP: DO J=OB%J1+1,OB%J2
            I_LOOP: DO I=OB%I1+1,OB%I2
               IC = M%CELL_INDEX(I,J,K)
               IF (.NOT.M%SOLID(IC)) CYCLE I_LOOP
               M%WALL_INDEX_HT3D(IC,:) = 0
               CELL_COUNT = 0

               MARCH_RIGHT: DO II=I+1,M%IBAR
                  ICN = M%CELL_INDEX(II,J,K)
                  IF (.NOT.M%SOLID(ICN)) THEN
                     CELL_COUNT(1) = II-I
                     M%WALL_INDEX_HT3D(IC,1) = M%WALL_INDEX(ICN,-1)
                     EXIT MARCH_RIGHT
                  ENDIF
               ENDDO MARCH_RIGHT

               MARCH_LEFT: DO II=I-1,1,-1
                  ICN = M%CELL_INDEX(II,J,K)
                  IF (.NOT.M%SOLID(ICN)) THEN
                     CELL_COUNT(-1) = I-II
                     M%WALL_INDEX_HT3D(IC,-1) = M%WALL_INDEX(ICN,1)
                     EXIT MARCH_LEFT
                  ENDIF
               ENDDO MARCH_LEFT

               MARCH_FORWARD: DO JJ=J+1,M%JBAR
                  ICN = M%CELL_INDEX(I,JJ,K)
                  IF (.NOT.M%SOLID(ICN)) THEN
                     CELL_COUNT(2) = JJ-J
                     M%WALL_INDEX_HT3D(IC,2) = M%WALL_INDEX(ICN,-2)
                     EXIT MARCH_FORWARD
                  ENDIF
               ENDDO MARCH_FORWARD

               MARCH_BACK: DO JJ=J-1,1,-1
                  ICN = M%CELL_INDEX(I,JJ,K)
                  IF (.NOT.M%SOLID(ICN)) THEN
                     CELL_COUNT(-2) = J-JJ
                     M%WALL_INDEX_HT3D(IC,-2) = M%WALL_INDEX(ICN,2)
                     EXIT MARCH_BACK
                  ENDIF
               ENDDO MARCH_BACK

               MARCH_UP: DO KK=K+1,M%KBAR
                  ICN = M%CELL_INDEX(I,J,KK)
                  IF (.NOT.M%SOLID(ICN)) THEN
                     CELL_COUNT(3) = KK-K
                     M%WALL_INDEX_HT3D(IC,3) = M%WALL_INDEX(ICN,-3)
                     EXIT MARCH_UP
                  ENDIF
               ENDDO MARCH_UP

               MARCH_DOWN: DO KK=K-1,1,-1
                  ICN = M%CELL_INDEX(I,J,KK)
                  IF (.NOT.M%SOLID(ICN)) THEN
                     CELL_COUNT(-3) = K-KK
                     M%WALL_INDEX_HT3D(IC,-3) = M%WALL_INDEX(ICN,3)
                     EXIT MARCH_DOWN
                  ENDIF
               ENDDO MARCH_DOWN

               ! Note: If multiple elements in the CELL_COUNT array have the same value, which will happen at a corner
               ! for example, then BACK=.TRUE. selects the last element.  This has the effect of giving precendence to +3
               ! at a (+1,+3) corner, for example.
               M%WALL_INDEX_HT3D(IC,0) = M%WALL_INDEX_HT3D(IC,MINLOC(CELL_COUNT,DIM=1,MASK=CELL_COUNT>0,BACK=.TRUE.)-4)

            ENDDO I_LOOP
         ENDDO J_LOOP
      ENDDO K_LOOP
   ENDIF HT3D_IF
ENDDO OBST_LOOP

END SUBROUTINE ASSIGN_HT3D_WALL_INDICES


!> \brief Determine if a point is in the interior of at least one mesh
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param ZZ z-coordinate of the point (m)

LOGICAL FUNCTION INTERIOR(XX,YY,ZZ)

INTEGER NM
REAL(EB), INTENT(IN) :: XX,YY,ZZ

INTERIOR = .FALSE.

DO NM=1,NMESHES
   IF (XX>MESHES(NM)%XS .AND. XX<MESHES(NM)%XF .AND. &
       YY>MESHES(NM)%YS .AND. YY<MESHES(NM)%YF .AND. &
       ZZ>MESHES(NM)%ZS .AND. ZZ<MESHES(NM)%ZF) INTERIOR = .TRUE.
ENDDO

END FUNCTION INTERIOR


!> \brief Assign the pressure zone for all cells connected to a given input point
!> \param NM Mesh number
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param ZZ z-coordinate of the point (m)
!> \param I_ZONE Index of the pressure zone to assign to the connected cells
!> \param I_ZONE_OVERLAP Index of a pressure zone that overlaps I_ZONE

SUBROUTINE ASSIGN_PRESSURE_ZONE(NM,XX,YY,ZZ,I_ZONE,I_ZONE_OVERLAP)

USE COMP_FUNCTIONS, ONLY : SHUTDOWN
REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(IN) :: NM,I_ZONE
INTEGER, INTENT(OUT) :: I_ZONE_OVERLAP
INTEGER :: NN,IOR,IC,IC_OLD,II,JJ,KK,III,JJJ,KKK,Q_N
INTEGER, ALLOCATABLE, DIMENSION(:) :: Q_I,Q_J,Q_K
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
TYPE (MESH_TYPE), POINTER :: M
TYPE (OBSTRUCTION_TYPE), POINTER :: OB

M=>MESHES(NM)
I_ZONE_OVERLAP = -1

IF (XX<M%XS-TWO_EPSILON_EB .OR. XX>M%XF+TWO_EPSILON_EB .OR. &
    YY<M%YS-TWO_EPSILON_EB .OR. YY>M%YF+TWO_EPSILON_EB .OR. &
    ZZ<M%ZS-TWO_EPSILON_EB .OR. ZZ>M%ZF+TWO_EPSILON_EB) RETURN

ALLOCATE(Q_I(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_J(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_K(M%IBAR*M%JBAR*M%KBAR))

! Find the cell indices corresponding to the given point

XI  = MAX( 1._EB , MIN( REAL(M%IBAR,EB)+ALMOST_ONE , M%CELLSI(NINT((XX-M%XS)*M%RDXINT)) + 1._EB ) )
YJ  = MAX( 1._EB , MIN( REAL(M%JBAR,EB)+ALMOST_ONE , M%CELLSJ(NINT((YY-M%YS)*M%RDYINT)) + 1._EB ) )
ZK  = MAX( 1._EB , MIN( REAL(M%KBAR,EB)+ALMOST_ONE , M%CELLSK(NINT((ZZ-M%ZS)*M%RDZINT)) + 1._EB ) )
II  = FLOOR(XI)
JJ  = FLOOR(YJ)
KK  = FLOOR(ZK)
IC  = M%CELL_INDEX(II,JJ,KK)

IF (M%SOLID(IC)) THEN
   WRITE(MESSAGE,'(A,I3,A)')  'ERROR: XYZ point for ZONE ',I_ZONE,' is inside a solid obstruction. Choose another XYZ.'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
   RETURN
ENDIF

! Add the first entry to "queue" of cells that need a pressure zone number

Q_I(1) = II
Q_J(1) = JJ
Q_K(1) = KK
Q_N    = 1
M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE

! Look to all cells adjacent to the starting cell and determine if they are in the ZONE as well.
! Repeat process until all cells are found.

SORT_QUEUE: DO

   IF (Q_N<1) EXIT SORT_QUEUE

   III = Q_I(Q_N)
   JJJ = Q_J(Q_N)
   KKK = Q_K(Q_N)
   Q_N = Q_N - 1

   SEARCH_LOOP: DO IOR=-3,3

      IF (IOR==0) CYCLE SEARCH_LOOP

      IC  = M%CELL_INDEX(III,JJJ,KKK)

      SELECT CASE(IOR)
         CASE(-1)
            IF (III==1) CYCLE SEARCH_LOOP
            II = III-1
            JJ = JJJ
            KK = KKK
         CASE( 1)
            IF (III==M%IBAR) CYCLE SEARCH_LOOP
            II = III+1
            JJ = JJJ
            KK = KKK
         CASE(-2)
            IF (JJJ==1) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ-1
            KK = KKK
         CASE( 2)
            IF (JJJ==M%JBAR) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ+1
            KK = KKK
         CASE(-3)
            IF (KKK==1) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ
            KK = KKK-1
         CASE( 3)
            IF (KKK==M%KBAR) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ
            KK = KKK+1
      END SELECT

      ! Look for thin obstructions bordering the current cell

      DO NN=1,M%N_OBST
         OB=>M%OBSTRUCTION(NN)
         SELECT CASE(IOR)
            CASE(-1)
               IF (II==  OB%I1 .AND. II==  OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE( 1)
               IF (II-1==OB%I1 .AND. II-1==OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE(-2)
               IF (JJ==  OB%J1 .AND. JJ==  OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE( 2)
               IF (JJ-1==OB%J1 .AND. JJ-1==OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE(-3)
               IF (KK==  OB%K1 .AND. KK==  OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
            CASE( 3)
               IF (KK-1==OB%K1 .AND. KK-1==OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
         END SELECT
      ENDDO

      ! If the last cell was solid and the current cell is not solid, stop the current directional march.

      IC_OLD = IC
      IC = M%CELL_INDEX(II,JJ,KK)

      IF (M%SOLID(IC_OLD) .AND. .NOT.M%SOLID(IC)) CYCLE SEARCH_LOOP

      IF (CC_IBM) THEN
         ! Here IBM_CGSC=1, IBM_SOLID=1:
         IF(M%CCVAR(III,JJJ,KKK,1)==1 .AND. M%CCVAR(II,JJ,KK,1)/=1) CYCLE SEARCH_LOOP

      ! If the current cell is not solid, but it is assigned another ZONE index, mark it as an overlap error and return

         ! Cell not SOLID for OBSTS, or GEOM cell not IBM_SOLID:
         IF (.NOT.M%SOLID(IC) .AND. M%CCVAR(II,JJ,KK,1)/=1 .AND. &
            M%PRESSURE_ZONE(II,JJ,KK)>=0 .AND.  M%PRESSURE_ZONE(II,JJ,KK)/=I_ZONE) THEN
            I_ZONE_OVERLAP = M%PRESSURE_ZONE(II,JJ,KK)
            RETURN
         ENDIF
      ELSE
         IF (.NOT.M%SOLID(IC) .AND. M%PRESSURE_ZONE(II,JJ,KK)>=0 .AND.  M%PRESSURE_ZONE(II,JJ,KK)/=I_ZONE) THEN
            I_ZONE_OVERLAP = M%PRESSURE_ZONE(II,JJ,KK)
            RETURN
         ENDIF
      ENDIF

      ! If the current cell is unassigned, assign the cell the ZONE index, I_ZONE, and then add this cell to the
      ! queue so that further searches might originate from it.

      IF (M%PRESSURE_ZONE(II,JJ,KK)<0) THEN
         M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
         Q_N      = Q_N+1
         Q_I(Q_N) = II
         Q_J(Q_N) = JJ
         Q_K(Q_N) = KK
      ENDIF

   ENDDO SEARCH_LOOP

ENDDO SORT_QUEUE

DEALLOCATE(Q_I)
DEALLOCATE(Q_J)
DEALLOCATE(Q_K)

END SUBROUTINE ASSIGN_PRESSURE_ZONE


!> \brief Determine the number of 1-D cells in a layer of solid material
!> \param DIFFUSIVITY \f$ k/(\rho c) \; (\hbox{m}^2/\hbox{s}) \f$, used to determine cell size
!> \param LAYER_THICKNESS Thickness of the material layer (m)
!> \param STRETCH_FACTOR A real number indicating the amount to stretch the second and subsequent cells in the interior
!> \param CELL_SIZE_FACTOR A real number indicating how much to shrink the first cell at the surface
!> \param N_LAYER_CELLS_MAX Maximum number of cells to assign to the layer
!> \param N_CELLS Number of cells in the layer
!> \param DX_MIN Minimum cell size

SUBROUTINE GET_N_LAYER_CELLS(DIFFUSIVITY,LAYER_THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR,N_LAYER_CELLS_MAX,N_CELLS,DX_MIN)

INTEGER, INTENT(OUT)  :: N_CELLS
INTEGER, INTENT(IN) :: N_LAYER_CELLS_MAX
REAL(EB), INTENT(OUT) :: DX_MIN
REAL(EB), INTENT(IN)  :: DIFFUSIVITY,LAYER_THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR
REAL(EB) :: DDSUM
INTEGER  :: N, I

DX_MIN = 0._EB
IF (LAYER_THICKNESS<=TWO_EPSILON_EB) THEN
   N_CELLS = 0
   DX_MIN = 0._EB
   RETURN
ENDIF

SHRINK_LOOP: DO N=1,N_LAYER_CELLS_MAX-1
   DDSUM = 0._EB
   SUM_LOOP: DO I=1,N
      DDSUM = DDSUM + STRETCH_FACTOR**(MIN(I-1,N-I))
   ENDDO SUM_LOOP
   IF ((LAYER_THICKNESS/DDSUM < CELL_SIZE_FACTOR*SQRT(DIFFUSIVITY)) .OR. (N==N_LAYER_CELLS_MAX-1)) THEN
      N_CELLS = N
      DX_MIN = LAYER_THICKNESS/DDSUM
      EXIT SHRINK_LOOP
   ENDIF
ENDDO SHRINK_LOOP

END SUBROUTINE GET_N_LAYER_CELLS


!> \brief Determine the coordinates of the 1-D interior cells
!> \param N_CELLS Number of cells in the entire collection of layers
!> \param N_LAYERS Number of layers
!> \param N_LAYER_CELLS Array holding the number of cells in each layer
!> \param SMALLEST_CELL_SIZE Array holding the sizes of the bounding cells in each layer (m)
!> \param STRETCH_FACTOR Array holding the stretching factors of each layer
!> \param X_S Array containing boundaries of the interior cells

SUBROUTINE GET_WALL_NODE_COORDINATES(N_CELLS,N_LAYERS,N_LAYER_CELLS,SMALLEST_CELL_SIZE,STRETCH_FACTOR,X_S)

INTEGER, INTENT(IN)     :: N_CELLS,N_LAYERS, N_LAYER_CELLS(N_LAYERS)
REAL(EB), INTENT(IN)    :: SMALLEST_CELL_SIZE(N_LAYERS),STRETCH_FACTOR(N_LAYERS)
REAL(EB), INTENT(OUT)   :: X_S(0:N_CELLS)

INTEGER I, II, NL
REAL(EB) DX_S

II = 0
X_S(0) = 0._EB
DO NL=1,N_LAYERS
   DO I=1,N_LAYER_CELLS(NL)
      II = II+1
      DX_S = SMALLEST_CELL_SIZE(NL)*STRETCH_FACTOR(NL)**(MIN(I-1,N_LAYER_CELLS(NL)-I))
      X_S(II) = X_S(II-1) + DX_S
   ENDDO
ENDDO

END SUBROUTINE GET_WALL_NODE_COORDINATES



!> \brief Determine the internal coordinates of the 1-D grid inside a solid.
!> \param N_CELLS Number of cells in the entire 1-D array
!> \param N_LAYERS Number of layers
!> \param N_LAYER_CELLS Array holding the number of cells in each layer
!> \param LAYER_THICKNESS Array of layer thicknesses (m)
!> \param GEOMETRY Indicates whether the surface is Cartesian, cylindrical, or spherical
!> \param X_S Array holding node coordinates (m)
!> \param LAYER_DIVIDE Number of layers that off gas to the front surface
!> \param DX Array of cell sizes (m)
!> \param RDX Reciprocal of DX (1/m)
!> \param RDXN Reciprocal of the distance from cell center to cell center (1/m)
!> \param DX_WGT Ratio of cell sizes: DX_WGT(I)=DX(I)*RDXN(I)*2
!> \param DXF Size of cell nearest the front surface (m)
!> \param DXB Size of cell nearest the back surface (m)
!> \param LAYER_INDEX Array of indices indicating the layer to which each interior cell belongs
!> \param MF_FRAC Array containing the fraction of each cells mass that is assigned to the front surface
!> \param INNER_RADIUS Inner radius of hollow cylinder or sphere (m)

SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
         LAYER_THICKNESS,GEOMETRY,X_S,LAYER_DIVIDE,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX,MF_FRAC,INNER_RADIUS)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS, N_LAYERS, N_LAYER_CELLS(N_LAYERS),GEOMETRY
REAL(EB), INTENT(IN)    :: X_S(0:N_CELLS),LAYER_THICKNESS(1:N_LAYERS),LAYER_DIVIDE,INNER_RADIUS
INTEGER, INTENT(OUT)    :: LAYER_INDEX(0:N_CELLS+1)
REAL(EB), INTENT(OUT)   :: DX(1:N_CELLS),RDX(0:N_CELLS+1),RDXN(0:N_CELLS),DX_WGT(0:N_CELLS),DXF,DXB, &
                           MF_FRAC(1:N_CELLS)

INTEGER :: I, II, NL, I_GRAD=0
REAL(EB) :: R, THICKNESS, X_DIVIDE

   THICKNESS = SUM(LAYER_THICKNESS)

   SELECT CASE(GEOMETRY)
      CASE(SURF_CARTESIAN)
         I_GRAD = 0
      CASE(SURF_CYLINDRICAL)
         I_GRAD = 1
      CASE(SURF_SPHERICAL)
         I_GRAD = 2
   END SELECT

   II = 0
   DO NL=1,N_LAYERS
      DO I=1,N_LAYER_CELLS(NL)
         II = II + 1
         LAYER_INDEX(II) = NL
      ENDDO
   ENDDO
   LAYER_INDEX(0) = 1
   LAYER_INDEX(N_CELLS+1) = N_LAYERS
   DXF = X_S(1)       - X_S(0)
   DXB = X_S(N_CELLS) - X_S(N_CELLS-1)

   ! Compute dx_weight for each node (dx_weight is the ratio of cell size to the combined size of the current and next cell)

   DO I=1,N_CELLS-1
      DX_WGT(I) = (X_S(I)-X_S(I-1))/(X_S(I+1)-X_S(I-1))
   ENDDO
   DX_WGT(0)       = 0.5_EB
   DX_WGT(N_CELLS) = 0.5_EB

   ! Compute dx and 1/dx for each node (dx is the distance from cell boundary to cell boundary)

   DO I=1,N_CELLS
      DX(I)  = X_S(I)-X_S(I-1)
      RDX(I) = 1._EB/DX(I)
   ENDDO

   ! Adjust 1/dx_n to 1/rdr for cylindrical case and 1/r^2dr for spaherical

   SELECT CASE(GEOMETRY)
      CASE(SURF_CYLINDRICAL)
         DO I=1,N_CELLS
            R = 0.5_EB*RDX(I)*((INNER_RADIUS+THICKNESS-X_S(I-1))**2-(INNER_RADIUS+THICKNESS-X_S(I))**2)
            RDX(I) = RDX(I)/R**I_GRAD
         ENDDO
      CASE(SURF_SPHERICAL)
         DO I=1,N_CELLS
            R = SQRT(ONTH*RDX(I)*((INNER_RADIUS+THICKNESS-X_S(I-1))**3-(INNER_RADIUS+THICKNESS-X_S(I))**3))
            RDX(I) = RDX(I)/R**I_GRAD
         ENDDO
   END SELECT
   RDX(0)         = RDX(1)
   RDX(N_CELLS+1) = RDX(N_CELLS)

   ! Compute 1/dx_n for each node (dx_n is the distance from cell center to cell center)

   DO I=1,N_CELLS-1
      RDXN(I) = 2._EB/(X_S(I+1)-X_S(I-1))
   ENDDO
   RDXN(0)       = 1._EB/(X_S(1)-X_S(0))
   RDXN(N_CELLS) = 1._EB/(X_S(N_CELLS)-X_S(N_CELLS-1))

   ! Adjust 1/dx_n to r/dr for cylindrical case and r^2/dr for spaherical

   IF (GEOMETRY /= SURF_CARTESIAN) THEN
      DO I=0,N_CELLS
         R = INNER_RADIUS+THICKNESS-X_S(I)
         RDXN(I) = RDXN(I)*R**I_GRAD
      ENDDO
   ENDIF

   ! Compute mass flux fraction array (array numbers indicate the fraction of mass flux that is added to the front

   IF (LAYER_DIVIDE >= REAL(N_LAYERS,EB)) THEN
      MF_FRAC = 1.0_EB
   ELSE
      MF_FRAC = 0._EB

      X_DIVIDE = 0._EB
      DO NL=1,N_LAYERS
         IF (LAYER_DIVIDE>=REAL(NL,EB)) THEN
            X_DIVIDE  = X_DIVIDE + LAYER_THICKNESS(NL)
         ELSE
            X_DIVIDE  = X_DIVIDE + MOD(LAYER_DIVIDE,1.0_EB)*LAYER_THICKNESS(NL)
            EXIT
         ENDIF
      ENDDO

      II = 0
      DIVILOOP: DO NL=1,N_LAYERS
         DO I=1,N_LAYER_CELLS(NL)
            II = II + 1
            IF (X_S(II) < X_DIVIDE) THEN
               MF_FRAC(II) = 1._EB
            ELSEIF (X_S(II-1) < X_DIVIDE) THEN
               MF_FRAC(II) = (X_DIVIDE-X_S(II-1))/DX(II)
               EXIT DIVILOOP
            ENDIF
         ENDDO
      ENDDO DIVILOOP
   ENDIF

END SUBROUTINE GET_WALL_NODE_WEIGHTS


!> \brief Determine weighting factors for 1-D solid cells
!> \param GEOMETRY Indicator of surface geometry: Cartesian, cylindrical, or spherical
!> \param NWP Number of interior cells
!> \param NWP_NEW Number of interior cells after shrinkage or swelling
!> \param INNER_RADIUS Inner radius of hollow cylinder or sphere (m)
!> \param X_S Array of interior cell edge positions (m)
!> \param X_S_NEW Array of interior cell edge positions after shrinkage or swelling (m)
!> \param INT_WGT Array of weighting factors for new arrangement of interior cells

SUBROUTINE GET_INTERPOLATION_WEIGHTS(GEOMETRY,NWP,NWP_NEW,INNER_RADIUS,X_S,X_S_NEW,INT_WGT)

INTEGER, INTENT(IN)  :: GEOMETRY,NWP,NWP_NEW
REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW), INNER_RADIUS
REAL(EB), INTENT(OUT) :: INT_WGT(NWP_NEW,NWP)

REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,VOL_NEW,VOL,THICKNESS
INTEGER I_NEW, I_OLD, I_GRAD


SELECT CASE(GEOMETRY)
   CASE(SURF_CARTESIAN)
      I_GRAD = 1
   CASE(SURF_CYLINDRICAL)
      I_GRAD = 2
   CASE(SURF_SPHERICAL)
      I_GRAD = 3
END SELECT

I_OLD = 1
INT_WGT = 0._EB
THICKNESS = X_S(NWP)

POINT_LOOP: DO I_NEW=1,NWP_NEW
XLOW_NEW = X_S_NEW(I_NEW-1)
XUP_NEW = X_S_NEW(I_NEW)

   OLD_POINT_LOOP: DO
      XLOW = X_S(I_OLD-1)
      XUP = X_S(I_OLD)
      VOL = (THICKNESS+INNER_RADIUS-XLOW)**I_GRAD-(THICKNESS+INNER_RADIUS-XUP)**I_GRAD
      VOL_NEW = (THICKNESS+INNER_RADIUS-MAX(XLOW_NEW,XLOW))**I_GRAD-(THICKNESS+INNER_RADIUS-MIN(XUP_NEW,XUP))**I_GRAD
      IF (VOL > 0._EB) INT_WGT(I_NEW,I_OLD) = MAX(0._EB,MIN(1._EB,VOL_NEW/VOL))
      IF (XUP >= XUP_NEW .OR. I_OLD==NWP) EXIT OLD_POINT_LOOP
      I_OLD = I_OLD+1
   ENDDO OLD_POINT_LOOP

ENDDO POINT_LOOP

END SUBROUTINE GET_INTERPOLATION_WEIGHTS


!> \brief Interpolates old array of wall properties to reflect new wall noding
!>
!> \param N_CELLS Maximum of the previous or current number of wall nodes
!> \param NWP Previous number of wall nodes
!> \param NWP_NEW Current number of wall nodes
!> \param INT_WGT array of weighting factors mapping old wall nodes to new wall nodes
!> \param ARR On input this is the old array and on output this is the new array

SUBROUTINE INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,ARR)

INTEGER, INTENT(IN)  :: N_CELLS,NWP,NWP_NEW
REAL(EB), INTENT(IN) :: INT_WGT(NWP_NEW,NWP)
REAL(EB), INTENT(INOUT ):: ARR(N_CELLS)
REAL(EB) :: TMP(N_CELLS)

INTEGER I,J

TMP = ARR
ARR = 0._EB
DO I = 1,NWP_NEW
  DO J = 1,NWP
        ARR(I) = ARR(I) + INT_WGT(I,J) * TMP(J)
  ENDDO
ENDDO

END SUBROUTINE INTERPOLATE_WALL_ARRAY


!> \brief Choose a random point from within a rectangular volume
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param ZZ z-coordinate of the point (m)
!> \param X1 Lower x-coordinate of the volume (m)
!> \param X2 Upper x-coordinate of the volume (m)
!> \param Y1 Lower y-coordinate of the volume (m)
!> \param Y2 Upper y-coordinate of the volume (m)
!> \param Z1 Lower z-coordinate of the volume (m)
!> \param Z2 Upper z-coordinate of the volume (m)

SUBROUTINE RANDOM_RECTANGLE(XX,YY,ZZ,X1,X2,Y1,Y2,Z1,Z2)

REAL(EB), INTENT(IN) :: X1,X2,Y1,Y2,Z1,Z2
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: RN

CALL RANDOM_NUMBER(RN)
XX = X1 + RN*(X2-X1)
CALL RANDOM_NUMBER(RN)
YY = Y1 + RN*(Y2-Y1)
CALL RANDOM_NUMBER(RN)
ZZ = Z1 + RN*(Z2-Z1)

END SUBROUTINE RANDOM_RECTANGLE


!> \brief Choose a random point from within a vertically oriented cone
!> \param NM Mesh number
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param ZZ z-coordinate of the point (m)
!> \param X0 x-coordinate of the center of the code base (m)
!> \param Y0 y-coordinate of the center of the code base (m)
!> \param Z0 z-coordinate of the center of the code base (m)
!> \param RR0 Radius of the base (m)
!> \param HH0 Height of the cone (m)

SUBROUTINE RANDOM_CONE(NM,XX,YY,ZZ,X0,Y0,Z0,RR0,HH0)

INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: X0,Y0,Z0,RR0,HH0
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: THETA,RR,RN
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
SEARCH_LOOP: DO
   CALL RANDOM_NUMBER(RN)
   ZZ = Z0 + HH0*RN
   CALL RANDOM_NUMBER(RN)
   THETA = TWOPI*RN
   CALL RANDOM_NUMBER(RN)
   RR = SQRT(RN)*RR0
   XX = X0 + RR*COS(THETA)
   YY = Y0 + RR*SIN(THETA)
   IF (RR**2>(RR0*(1._EB-(ZZ-Z0)/HH0))**2) CYCLE SEARCH_LOOP
   IF (XX>=M%XS .AND. XX<=M%XF .AND. YY>=M%YS .AND. YY<=M%YF .AND. ZZ>=M%ZS .AND. ZZ<=M%ZF) EXIT SEARCH_LOOP
ENDDO SEARCH_LOOP

END SUBROUTINE RANDOM_CONE


!> \brief Choose a random point on a horizontally-oriented ring
!> \param NM Mesh number
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param X0 x-coordinate of the center of the ring (m)
!> \param Y0 y-coordinate of the center of the ring (m)
!> \param RR0 Radius of the ring (m)

SUBROUTINE RANDOM_RING(NM,XX,YY,X0,Y0,RR0)

INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: X0,Y0,RR0
REAL(EB), INTENT(OUT) :: XX,YY
REAL(EB) :: THETA,RN
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
SEARCH_LOOP: DO
   CALL RANDOM_NUMBER(RN)
   THETA = TWOPI*RN
   XX = X0 + RR0*COS(THETA)
   YY = Y0 + RR0*SIN(THETA)
   IF (XX>=M%XS .AND. XX<=M%XF .AND. YY>=M%YS .AND. YY<=M%YF) EXIT SEARCH_LOOP
ENDDO SEARCH_LOOP

END SUBROUTINE RANDOM_RING


!> \brief Assign an ordered point on a horizontally-oriented ring
!> \param NM Mesh number
!> \param XX x-coordinate of the point (m)
!> \param YY y-coordinate of the point (m)
!> \param X0 x-coordinate of the center of the ring (m)
!> \param Y0 y-coordinate of the center of the ring (m)
!> \param RR0 Radius of the ring (m)
!> \param NP Number of the current point
!> \param NP_TOTAL Total number of points in the ring

SUBROUTINE UNIFORM_RING(NM,XX,YY,X0,Y0,RR0,NP,NP_TOTAL)

INTEGER, INTENT(IN) :: NM,NP,NP_TOTAL
REAL(EB), INTENT(IN) :: X0,Y0,RR0
REAL(EB), INTENT(OUT) :: XX,YY
REAL(EB) :: THETA
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
SEARCH_LOOP: DO
   THETA = TWOPI*REAL(NP,EB)/REAL(NP_TOTAL,EB)
   XX = X0 + RR0*COS(THETA)
   YY = Y0 + RR0*SIN(THETA)
   IF (XX>=M%XS .AND. XX<=M%XF .AND. YY>=M%YS .AND. YY<=M%YF) EXIT SEARCH_LOOP
ENDDO SEARCH_LOOP

END SUBROUTINE UNIFORM_RING


!> \brief Estimate area of intersection of circle and rectangle
!> \param X0 x-coordinate of the center of the circle (m)
!> \param Y0 y-coordinate of the center of the circle (m)
!> \param RAD Radius of the circle (m)
!> \param X1 Lower x-coordinate of the rectangle (m)
!> \param X2 Upper x-coordinate of the rectangle (m)
!> \param Y1 Lower y-coordinate of the rectangle (m)
!> \param Y2 Upper y-coordinate of the rectangle (m)

REAL(EB) FUNCTION CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RAD,X1,X2,Y1,Y2)

REAL(EB), INTENT(IN) :: X0,Y0,RAD,X1,X2,Y1,Y2
INTEGER :: NN,II,JJ
REAL(EB) :: DELTA_AREA,XX,YY

CIRCLE_CELL_INTERSECTION_AREA = 0._EB
NN = 50
DELTA_AREA = (X2-X1)*(Y2-Y1)/REAL(NN,EB)**2

DO JJ=1,NN
   DO II=1,NN
      XX = X1 + (X2-X1)*(II-0.5_EB)/REAL(NN,EB)
      YY = Y1 + (Y2-Y1)*(JJ-0.5_EB)/REAL(NN,EB)
      IF ( ((XX-X0)**2+(YY-Y0)**2)<RAD**2 ) CIRCLE_CELL_INTERSECTION_AREA = CIRCLE_CELL_INTERSECTION_AREA + DELTA_AREA
   ENDDO
ENDDO

END FUNCTION CIRCLE_CELL_INTERSECTION_AREA


!> \brief Calculate volume of the interaction of MESH NM with a cone
!> \param NM Mesh number
!> \param X0 x-coordinate of the center of the base of the cone (m)
!> \param Y0 y-coordinate of the center of the base of the cone (m)
!> \param Z0 z-coordinate of the center of the base of the cone (m)
!> \param RR0 Radius of the base of the cone (m)
!> \param HH0 Height of the cone (m)

REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,HH0)

INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: X0,Y0,Z0,RR0,HH0
INTEGER :: I,J,K,NX,NY,NZ
REAL(EB) :: XX,YY,ZZ,R2,X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX,DX,DY,DZ
TYPE (MESH_TYPE), POINTER :: M

CONE_MESH_INTERSECTION_VOLUME = 0._EB
M => MESHES(NM)

X_MIN = MAX(M%XS,X0-RR0)
X_MAX = MIN(M%XF,X0+RR0)
Y_MIN = MAX(M%YS,Y0-RR0)
Y_MAX = MIN(M%YF,Y0+RR0)
Z_MIN = MAX(M%ZS,Z0)
Z_MAX = MIN(M%ZF,Z0+HH0)
IF (X_MAX<=X_MIN .OR. Y_MAX<=Y_MIN .OR. Z_MAX<=Z_MIN) RETURN
NX = CEILING(10*(X_MAX-X_MIN))
NY = CEILING(10*(Y_MAX-Y_MIN))
NZ = CEILING(10*(Z_MAX-Z_MIN))
DX = (X_MAX-X_MIN)/REAL(NX,EB)
DY = (Y_MAX-Y_MIN)/REAL(NY,EB)
DZ = (Z_MAX-Z_MIN)/REAL(NZ,EB)

DO K=1,NZ
   ZZ = Z_MIN + (K-0.5_EB)*DZ
   IF (ZZ<Z0 .OR. ZZ>Z0+HH0) CYCLE
   DO J=1,NY
      YY = Y_MIN + (J-0.5_EB)*DY
      DO I=1,NX
         XX = X_MIN + (I-0.5_EB)*DX
         R2 = (XX-X0)**2+(YY-Y0)**2
         IF (R2<(RR0*(1._EB-(ZZ-Z0)/HH0))**2) &
            CONE_MESH_INTERSECTION_VOLUME = CONE_MESH_INTERSECTION_VOLUME + DX*DY*DZ
      ENDDO
   ENDDO
ENDDO

END FUNCTION CONE_MESH_INTERSECTION_VOLUME


!> \brief Rotate and translate a point using the parameters of MOVE input line
!> \param X x-coordinate of the point (m)
!> \param Y y-coordinate of the point (m)
!> \param Z z-coordinate of the point (m)
!> \param MOVE_INDEX Index of the set of parameters on a MOVE input line

SUBROUTINE TRANSFORM_COORDINATES(X,Y,Z,MOVE_INDEX)

INTEGER, INTENT(IN) :: MOVE_INDEX
REAL(EB), INTENT(INOUT) :: X,Y,Z
REAL(EB) :: M(3,3),UP(3,1),S(3,3),UUT(3,3),IDENTITY(3,3),X_VECTOR(3,1),X_VECTOR_0(3,1),X_VECTOR_1(4,1),SCL(3,3)
TYPE(MOVEMENT_TYPE), POINTER :: MV

MV => MOVEMENT(MOVE_INDEX)

! If T34 present, use Matrix T3x4 for transformation: Note it takes precedence over other transformation parameters.

IF (SUM(ABS(MV%T34(1:3,1:4))) > TWO_EPSILON_EB) THEN
   X_VECTOR_1 = RESHAPE( (/X,Y,Z,1._EB/), (/4, 1/) )
   X_VECTOR   = MATMUL(MV%T34, X_VECTOR_1)
   X = X_VECTOR(1,1)
   Y = X_VECTOR(2,1)
   Z = X_VECTOR(3,1)
   RETURN
ENDIF

X_VECTOR = RESHAPE( (/X,Y,Z/),(/3,1/) )
X_VECTOR_0 = RESHAPE( (/MV%X0,MV%Y0,MV%Z0/),(/3,1/) )
UP = RESHAPE(MV%AXIS,(/3,1/))
S =  RESHAPE( (/  0.0_EB, -UP(3,1),  UP(2,1),&
                 UP(3,1),   0.0_EB, -UP(1,1),&
                -UP(2,1),  UP(1,1),  0.0_EB  /),(/3,3/))
UUT = MATMUL(UP,TRANSPOSE(UP))
IDENTITY = RESHAPE ((/ 1.0_EB,0.0_EB,0.0_EB,&
                       0.0_EB,1.0_EB,0.0_EB,&
                       0.0_EB,0.0_EB,1.0_EB /),(/3,3/))
M = UUT + COS(MV%ROTATION_ANGLE*DEG2RAD)*(IDENTITY - UUT) + SIN(MV%ROTATION_ANGLE*DEG2RAD)*S

SCL= RESHAPE((/ MV%SCALE*MV%SCALEX,0.0_EB,0.0_EB, &
                0.0_EB,MV%SCALE*MV%SCALEY,0.0_EB, &
                0.0_EB,0.0_EB,MV%SCALE*MV%SCALEZ /),(/3,3/))

! Scaling takes precedence over rotation transformation (applied first on the local axes):

X_VECTOR = X_VECTOR_0 + MATMUL(MATMUL(M,SCL),X_VECTOR-X_VECTOR_0)
X = X_VECTOR(1,1) + MV%DX
Y = X_VECTOR(2,1) + MV%DY
Z = X_VECTOR(3,1) + MV%DZ

END SUBROUTINE TRANSFORM_COORDINATES

END MODULE GEOMETRY_FUNCTIONS


!> \brief Routines that do various mathematical manipulations

MODULE MATH_FUNCTIONS

USE PRECISION_PARAMETERS
IMPLICIT NONE

CONTAINS


!> \brief Add new value to a growing histogram
!> \param NBINS Number of bins in the histogram
!> \param LIMITS The range of values underlying the histogram
!> \param COUNTS Array of weight currently assigned to each bin
!> \param VAL New value to be added
!> \param WEIGHT Relative weight associated with the new value

SUBROUTINE UPDATE_HISTOGRAM(NBINS,LIMITS,COUNTS,VAL,WEIGHT)
INTEGER,INTENT(IN)::NBINS
REAL(EB), INTENT(IN)::LIMITS(2),VAL,WEIGHT
REAL(EB), INTENT(INOUT) :: COUNTS(NBINS)
INTEGER::IND=0
IND=MIN(NBINS,MAX(CEILING((VAL-LIMITS(1))/(LIMITS(2)-LIMITS(1))*NBINS),1))
COUNTS(IND)=COUNTS(IND)+WEIGHT
END SUBROUTINE UPDATE_HISTOGRAM


!> \brief Linearly interpolate the value at a point in a 3D array
!> \param A The 3D array of values
!> \param I The lower x index of the array
!> \param J The lower y index of the array
!> \param K The lower z index of the array
!> \param P Fraction of the distance from the lower to upper x coordinate
!> \param R Fraction of the distance from the lower to upper y coordinate
!> \param S Fraction of the distance from the lower to upper z coordinate

REAL(EB) FUNCTION AFILL2(A,I,J,K,P,R,S)

REAL(EB), INTENT(IN), DIMENSION(0:,0:,0:) :: A
INTEGER, INTENT(IN) :: I,J,K
REAL(EB) A111,A211,A121,A221,A112,A212,A122,A222,P,R,S,PP,RR,SS

A111 = A(I,J,K)
A211 = A(I+1,J,K)
A121 = A(I,J+1,K)
A221 = A(I+1,J+1,K)
A112 = A(I,J,K+1)
A212 = A(I+1,J,K+1)
A122 = A(I,J+1,K+1)
A222 = A(I+1,J+1,K+1)
PP = 1._EB-P
RR = 1._EB-R
SS = 1._EB-S
AFILL2 = ((PP*A111+P*A211)*RR+(PP*A121+P*A221)*R)*SS+ ((PP*A112+P*A212)*RR+(PP*A122+P*A222)*R)*S

END FUNCTION AFILL2


!> \brief Calculate the value of polynomial function.
!> \param N Number of coefficients in the polynomial
!> \param TEMP The independent variable
!> \param COEF The array of coefficients

REAL(EB) FUNCTION POLYVAL(N,TEMP,COEF)
INTEGER, INTENT(IN) :: N
REAL(EB), INTENT(IN) :: TEMP,COEF(N)
INTEGER :: I
POLYVAL = 0._EB
DO I=1,N
   POLYVAL  = POLYVAL  + COEF(I)*TEMP**(I-1)
ENDDO
END FUNCTION POLYVAL


!> \brief Determine the index of a RAMP with a given name (ID)
!> \param ID The name of the ramp function
!> \param TYPE The kind of ramp
!> \param RAMP_INDEX The index of the ramp
!> \param DUPLICATE_RAMP Optional logical parameter indicating if the ramp has already been processed

SUBROUTINE GET_RAMP_INDEX(ID,TYPE,RAMP_INDEX,DUPLICATE_RAMP)

USE GLOBAL_CONSTANTS, ONLY: N_RAMP,RAMP_ID,RAMP_TYPE
USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_CHARACTER_ARRAY
CHARACTER(*), INTENT(IN) :: ID,TYPE
INTEGER, INTENT(OUT) :: RAMP_INDEX
INTEGER :: NR
LOGICAL, INTENT(IN), OPTIONAL :: DUPLICATE_RAMP

IF (ID=='null') THEN
   RAMP_INDEX = 0
   RETURN
ENDIF

IF (.NOT.PRESENT(DUPLICATE_RAMP)) THEN
   SEARCH: DO NR=1,N_RAMP
      IF (ID==RAMP_ID(NR)) THEN
         RAMP_INDEX = NR
         RETURN
      ENDIF
   ENDDO SEARCH
ENDIF

IF (N_RAMP>=SIZE(RAMP_ID)) THEN
   RAMP_ID   => REALLOCATE_CHARACTER_ARRAY(RAMP_ID,  LABEL_LENGTH,1,SIZE(RAMP_ID)+100)
   RAMP_TYPE => REALLOCATE_CHARACTER_ARRAY(RAMP_TYPE,LABEL_LENGTH,1,SIZE(RAMP_ID)+100)
ENDIF

N_RAMP                = N_RAMP + 1
RAMP_INDEX            = N_RAMP
RAMP_ID(RAMP_INDEX)   = ID
RAMP_TYPE(RAMP_INDEX) = TYPE

END SUBROUTINE GET_RAMP_INDEX


!> \brief Determine the index of a table with a given name
!> \param ID Name of the table
!> \param TYPE Kind of table
!> \param TABLE_INDEX Index of the table

SUBROUTINE GET_TABLE_INDEX(ID,TYPE,TABLE_INDEX)

USE GLOBAL_CONSTANTS, ONLY: N_TABLE,TABLE_ID,TABLE_TYPE
CHARACTER(*), INTENT(IN) :: ID
INTEGER, INTENT(IN) :: TYPE
INTEGER, INTENT(OUT) :: TABLE_INDEX
INTEGER :: NT

IF (ID=='null') THEN
   TABLE_INDEX = 0
   RETURN
ENDIF

SEARCH: DO NT=1,N_TABLE
   IF (ID==TABLE_ID(NT)) THEN
      TABLE_INDEX = NT
      RETURN
   ENDIF
ENDDO SEARCH

N_TABLE                = N_TABLE + 1
TABLE_INDEX            = N_TABLE
TABLE_ID(TABLE_INDEX)   = ID
TABLE_TYPE(TABLE_INDEX) = TYPE

END SUBROUTINE GET_TABLE_INDEX


!> \brief Return an interpolated value from a given RAMP function
!> \param RAMP_INPUT Independent variable of the RAMP function
!> \param TAU Time scale used for reserved t-squared or tanh ramps
!> \param RAMP_INDEX Index of the RAMP function

REAL(EB) FUNCTION EVALUATE_RAMP(RAMP_INPUT,TAU,RAMP_INDEX)

USE TYPES, ONLY: RAMPS
USE DEVICE_VARIABLES, ONLY: DEVICE
USE CONTROL_VARIABLES, ONLY: CONTROL
REAL(EB), INTENT(IN) :: RAMP_INPUT,TAU
REAL(EB):: RAMP_POSITION
INTEGER, INTENT(IN)  :: RAMP_INDEX

SELECT CASE(RAMP_INDEX)
   CASE(-2)
      EVALUATE_RAMP = MAX( TANH(RAMP_INPUT/TAU), 0._EB )
   CASE(-1)
      EVALUATE_RAMP = MIN( (RAMP_INPUT/TAU)**2 , 1.0_EB )
   CASE( 0)
      EVALUATE_RAMP = 1._EB
   CASE(1:)
      IF (RAMPS(RAMP_INDEX)%DEVC_INDEX > 0) THEN
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,DEVICE(RAMPS(RAMP_INDEX)%DEVC_INDEX)%SMOOTHED_VALUE - RAMPS(RAMP_INDEX)%T_MIN))
      ELSEIF (RAMPS(RAMP_INDEX)%CTRL_INDEX > 0) THEN
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,CONTROL(RAMPS(RAMP_INDEX)%CTRL_INDEX)%INSTANT_VALUE - RAMPS(RAMP_INDEX)%T_MIN))
      ELSE
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,RAMP_INPUT - RAMPS(RAMP_INDEX)%T_MIN))
      ENDIF
      EVALUATE_RAMP = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(NINT(RAMP_POSITION*RAMPS(RAMP_INDEX)%RDT))
END SELECT

END FUNCTION EVALUATE_RAMP


!> \brief Compute inverse erfc function
!> \param Y Y=ERFC(X), where X is returned by IERFC(Y)

REAL(EB) FUNCTION IERFC(Y)

REAL(EB), INTENT(IN) :: Y
REAL(EB) :: QA,QB,QC,QD,Q0,Q1,Q2,Q3,Q4,PA,PB,P0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10, &
            P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22,X,S,T,U,W,Z
PARAMETER ( &
QA = 9.16461398268964-01_EB, &
QB = 2.31729200323405-01_EB, &
QC = 4.88826640273108-01_EB, &
QD = 1.24610454613712-01_EB, &
Q0 = 4.99999303439796-01_EB, &
Q1 = 1.16065025341614-01_EB, &
Q2 = 1.50689047360223-01_EB, &
Q3 = 2.69999308670029-01_EB, &
Q4 = -7.28846765585675-02_EB)
PARAMETER ( &
PA = 3.97886080735226000+00_EB, &
PB = 1.20782237635245222-01_EB, &
P0 = 2.44044510593190935-01_EB, &
P1 = 4.34397492331430115-01_EB, &
P2 = 6.86265948274097816-01_EB, &
P3 = 9.56464974744799006-01_EB, &
P4 = 1.16374581931560831+00_EB, &
P5 = 1.21448730779995237+00_EB, &
P6 = 1.05375024970847138+00_EB, &
P7 = 7.13657635868730364-01_EB, &
P8 = 3.16847638520135944-01_EB, &
P9 = 1.47297938331485121-02_EB, &
P10 = -1.05872177941595488-01_EB, &
P11 = -7.43424357241784861-02_EB)
PARAMETER ( &
P12 = 2.20995927012179067-03_EB, &
P13 = 3.46494207789099922-02_EB, &
P14 = 1.42961988697898018-02_EB, &
P15 = -1.18598117047771104-02_EB, &
P16 = -1.12749169332504870-02_EB, &
P17 = 3.39721910367775861-03_EB, &
P18 = 6.85649426074558612-03_EB, &
P19 = -7.71708358954120939-04_EB, &
P20 = -3.51287146129100025-03_EB, &
P21 = 1.05739299623423047-04_EB, &
P22 = 1.12648096188977922-03_EB)

Z = Y
IF (Y  > 1._EB) Z = 2._EB - Y
W = QA - LOG(Z)
U = SQRT(W)
S = (QC + LOG(U)) / W
T = 1._EB / (U + QB)

X = U * (1._EB - S * (0.5_EB + S * QD)) - ((((Q4 * T + Q3) * T + Q2) * T + Q1) * T + Q0) * T
T = PA / (PA + X)
U = T - 0.5_EB

S = (((((((((P22 * U + P21) * U + P20) * U + P19) * U + P18) * U + P17) * U + P16) * U + P15) * U + P14) * U + P13) * U + P12

S = ((((((((((((S * U + P11) * U + P10) * U +  P9) * U + P8) * U + P7) * U + P6) * U + P5) * U + P4) * U + P3) &
    * U + P2) * U + P1) * U + P0) * T - Z * EXP(X * X - PB)

X = X + S * (1._EB + X * S)

IF (Y > 1._EB) X = -X

IERFC = X

END FUNCTION IERFC


!> \brief Solve a linear system of equations with Gauss-Jordon elimination (Press, Numerical Recipes)
!> \param A Primary matrix of A*x=b
!> \param N Dimension of A
!> \param NP Dimension of array containing A
!> \param B Array of vectors B of the linear system
!> \param M Number of columns of B
!> \param MP Number of columns of the array holding the B vectors
!> \param IERROR Error code

SUBROUTINE GAUSSJ(A,N,NP,B,M,MP,IERROR)

INTEGER, INTENT(IN) :: M,MP,N,NP
REAL(EB), INTENT(INOUT) :: A(NP,NP),B(NP,MP)
INTEGER, INTENT(OUT) :: IERROR
REAL(EB) :: BIG,DUM,PIVINV
INTEGER :: I,ICOL=0,IROW=0,J,K,L,LL,INDXC(NP),INDXR(NP),IPIV(NP)

IERROR = 0
IPIV(1:N) = 0

DO I=1,N
   BIG = 0._EB
   DO J=1,N
      IF (IPIV(J)/=1) THEN
         DO K=1,N
            IF (IPIV(K)==0) THEN
               IF (ABS(A(J,K))>=BIG) THEN
                  BIG = ABS(A(J,K))
                  IROW = J
                  ICOL = K
               ENDIF
            ELSE IF (IPIV(K)>1) THEN
               IERROR = 103   ! Singular matrix in gaussj
               RETURN
            ENDIF
         ENDDO
      ENDIF
   ENDDO
   IPIV(ICOL) = IPIV(ICOL) + 1
   IF (IROW/=ICOL) THEN
      DO L=1,N
         DUM = A(IROW,L)
         A(IROW,L) = A(ICOL,L)
         A(ICOL,L) = DUM
      ENDDO
      DO L=1,M
         DUM = B(IROW,L)
         B(IROW,L) = B(ICOL,L)
         B(ICOL,L) = DUM
      ENDDO
   ENDIF
   INDXR(I) = IROW
   INDXC(I) = ICOL
   IF (ABS(A(ICOL,ICOL))<=TWO_EPSILON_EB) THEN
      IERROR = 103  ! Singular matrix in gaussj
      RETURN
      ENDIF
   PIVINV = 1._EB/A(ICOL,ICOL)
   A(ICOL,ICOL) = 1._EB
   A(ICOL,1:N) = A(ICOL,1:N) * PIVINV
   B(ICOL,1:M) = B(ICOL,1:M) * PIVINV
   DO LL=1,N
      IF (LL/=ICOL) THEN
         DUM = A(LL,ICOL)
         A(LL,ICOL) = 0._EB
         A(LL,1:N) = A(LL,1:N) - A(ICOL,1:N)*DUM
         B(LL,1:M) = B(LL,1:M) - B(ICOL,1:M)*DUM
      ENDIF
   ENDDO
ENDDO
DO L=N,1,-1
   IF (INDXR(L)/=INDXC(L)) THEN
      DO K=1,N
         DUM = A(K,INDXR(L))
         A(K,INDXR(L)) = A(K,INDXC(L))
         A(K,INDXC(L)) = DUM
      ENDDO
   ENDIF
ENDDO

END SUBROUTINE GAUSSJ


!> \brief Linearly interpolate the value of a given function at a given point
!> \param X Independent variable
!> \param Y Dependent variable
!> \param XI Point to interpolate
!> \param ANS Value of the function at the point XI

SUBROUTINE INTERPOLATE1D(X,Y,XI,ANS)

REAL(EB), INTENT(IN), DIMENSION(:) :: X, Y
REAL(EB), INTENT(IN) :: XI
REAL(EB), INTENT(OUT) :: ANS
INTEGER I, UX,LX

UX = UBOUND(X,1)
LX = LBOUND(X,1)

IF (XI <= X(LX)) THEN
   ANS = Y(LX)
ELSEIF (XI >= X(UX)) THEN
   ANS = Y(UX)
ELSE
   L1: DO I=LX,UX-1
      IF (ABS(XI -X(I)) <= SPACING(X(I))) THEN
         ANS = Y(I)
         EXIT L1
      ELSEIF (X(I+1)>XI) THEN
         ANS = Y(I)+(XI-X(I))/(X(I+1)-X(I)) * (Y(I+1)-Y(I))
         EXIT L1
      ENDIF
   ENDDO L1
ENDIF

END SUBROUTINE INTERPOLATE1D


!> \brief Interpolate a 1D array of numbers
!> \param LOWER Lower index of the array X
!> \param X Array of numbers
!> \param XI Real number representing a fractional array index
!> \param ANS Interpolated value at XI

SUBROUTINE INTERPOLATE1D_UNIFORM(LOWER,X,XI,ANS)

INTEGER, INTENT(IN) :: LOWER
REAL(EB), INTENT(IN), DIMENSION(LOWER:) :: X
REAL(EB), INTENT(IN) :: XI
REAL(EB), INTENT(OUT) :: ANS
INTEGER I, UX,LX
REAL(EB) :: FRAC

UX = UBOUND(X,1)
LX = LBOUND(X,1)

IF (XI <= LX) THEN
   ANS = X(LX)
ELSEIF (XI >= UX) THEN
   ANS = X(UX)
ELSE
   I = INT(XI)
   FRAC = XI - REAL(I,EB)
   ANS = X(I) + FRAC*(X(I+1)-X(I))
ENDIF

END SUBROUTINE INTERPOLATE1D_UNIFORM


!> \brief Find interpolated value in a 2D array
!> \param TABLE_INDEX Index of the table
!> \param XI Index of first array dimension
!> \param YI Index of second array dimension
!> \param ANS Interpolated value at (XI,YI)

SUBROUTINE INTERPOLATE2D(TABLE_INDEX,XI,YI,ANS)

USE TYPES, ONLY: TABLES,TABLES_TYPE
INTEGER, INTENT(IN) :: TABLE_INDEX
REAL(EB), INTENT(IN) :: XI,YI
REAL(EB), INTENT(OUT) :: ANS
REAL(EB) :: XL,XU,FRAC
INTEGER :: I,J
TYPE (TABLES_TYPE), POINTER :: TA

TA => TABLES(TABLE_INDEX)

! Do 1D for X edges of table
IF (XI <= TA%LX) THEN
   CALL INTERPOLATE1D(TA%Y,TA%Z(1,:),YI,ANS)
   RETURN
ELSEIF (XI >= TA%UX) THEN
   CALL INTERPOLATE1D(TA%Y,TA%Z(TA%NUMBER_ROWS,:),YI,ANS)
   RETURN
ENDIF

! Do 1D for Y edges of table
IF (YI <= TA%LY) THEN
   CALL INTERPOLATE1D(TA%X,TA%Z(:,1),XI,ANS)
   RETURN
ELSEIF(YI >= TA%UY) THEN
   CALL INTERPOLATE1D(TA%X,TA%Z(:,TA%NUMBER_COLUMNS),XI,ANS)
   RETURN
ENDIF

! Search for column, do 1D if exact
DO J = 1, TA%NUMBER_COLUMNS
   IF (ABS(YI-TA%Y(J)) < SPACING(TA%Y(J))) THEN
      CALL INTERPOLATE1D(TA%X,TA%Z(:,J),XI,ANS)
      RETURN
   ENDIF
   IF (TA%Y(J) > YI) EXIT
ENDDO

! Search for row, do 1D if exact
DO I = 1, TA%NUMBER_ROWS
   IF(ABS(XI-TA%X(I)) < SPACING(TA%X(I))) THEN
      CALL INTERPOLATE1D(TA%Y,TA%Z(I,:),YI,ANS)
      RETURN
   ENDIF
   IF (TA%X(I) > XI) EXIT
ENDDO

FRAC = (XI-TA%X(I-1))/(TA%X(I)-TA%X(I-1))
XL = FRAC*(TA%Z(I,J-1)-TA%Z(I-1,J-1))+TA%Z(I-1,J-1)
XU = FRAC*(TA%Z(I,J  )-TA%Z(I-1,J  ))+TA%Z(I-1,J)
ANS = (YI-TA%Y(J-1))/(TA%Y(J)-TA%Y(J-1))*(XU-XL)+XL

END SUBROUTINE INTERPOLATE2D


!> \brief Randomly choose a point from a normal distribution
!> \param MEAN Mean of the normal distribution
!> \param SIGMA Standard deviation

REAL(EB) FUNCTION NORMAL(MEAN,SIGMA)

REAL(EB), INTENT(IN) :: MEAN,SIGMA
REAL(EB) :: TMP,FAC,GSAVE,RSQ,R1,R2
REAL     :: RN
INTEGER :: FLAG
SAVE FLAG,GSAVE
DATA FLAG /0/

IF (FLAG==0) THEN
   RSQ=2.0_EB
   DO WHILE(RSQ>=1.0_EB.OR.RSQ==0.0_EB)
      CALL RANDOM_NUMBER(RN)
      R1=2.0_EB*REAL(RN,EB)-1.0_EB
      CALL RANDOM_NUMBER(RN)
      R2=2.0_EB*REAL(RN,EB)-1.0_EB
      RSQ=R1*R1+R2*R2
   ENDDO
   FAC=SQRT(-2.0_EB*LOG(RSQ)/RSQ)
   GSAVE=R1*FAC
   TMP=R2*FAC
   FLAG=1
ELSE
   TMP=GSAVE
   FLAG=0
ENDIF
NORMAL=TMP*SIGMA+MEAN

END FUNCTION NORMAL


!> \brief Compute the cross product of two triplets, A x B = C
!> \param C The resulting vector
!> \param A First vector
!> \param B Second vector

SUBROUTINE CROSS_PRODUCT(C,A,B)

REAL(EB), INTENT(IN) :: A(3),B(3)
REAL(EB), INTENT(OUT) :: C(3)

C(1) = A(2)*B(3)-A(3)*B(2)
C(2) = A(3)*B(1)-A(1)*B(3)
C(3) = A(1)*B(2)-A(2)*B(1)

END SUBROUTINE CROSS_PRODUCT


!> \brief Randomly choose a value from the distribution with given CDF
!> \param CDF Cumulative Distribution Function
!> \param VAR Independent variable
!> \param NPTS Number of points in the CDF
!> \param CHOICE Randomly chosen value

SUBROUTINE RANDOM_CHOICE(CDF,VAR,NPTS,CHOICE)

INTEGER,  INTENT(IN)  :: NPTS
REAL(EB), INTENT(IN)  :: CDF(0:NPTS),VAR(0:NPTS)
REAL(EB), INTENT(OUT) :: CHOICE
INTEGER  :: IT
REAL(EB) :: CFRAC,A,B
REAL(EB) :: RN
REAL     :: RN2

CALL RANDOM_NUMBER(RN2)
RN = REAL(RN2,EB)
A = MINVAL(CDF)
B = MAXVAL(CDF)
RN = A + (B-A)*RN

CDF_LOOP: DO IT=1,NPTS
   IF (CDF(IT) > RN) THEN
      CFRAC  = (RN-CDF(IT-1))/(CDF(IT)-CDF(IT-1))
      CHOICE = VAR(IT-1) + (VAR(IT)-VAR(IT-1))*CFRAC
      EXIT CDF_LOOP
   ENDIF
ENDDO CDF_LOOP

END SUBROUTINE RANDOM_CHOICE


REAL(EB) FUNCTION MINMOD2(X,Y)
REAL(EB), INTENT(IN) :: X,Y
MINMOD2 = 0.5_EB*(SIGN(1._EB,X)+SIGN(1._EB,Y))*MIN(ABS(X),ABS(Y))
END FUNCTION MINMOD2


REAL(EB) FUNCTION MINMOD4(W,X,Y,Z)
REAL(EB), INTENT(IN) :: W,X,Y,Z
MINMOD4 = 0.125_EB*(SIGN(1._EB,W)+SIGN(1._EB,X))* &
          ABS( (SIGN(1._EB,W)+SIGN(1._EB,Y))*(SIGN(1._EB,W)+SIGN(1._EB,Z)) )* &
          MIN(ABS(W),ABS(X),ABS(Y),ABS(Z))
END FUNCTION MINMOD4


!> \brief Generate pairs of normally distributed pseudo-random numbers with zero mean and unit variance based on the
!> Box-Muller transformation.
!> \param Z0 Output value 1
!> \param Z1 Output value 2

SUBROUTINE BOX_MULLER(Z0,Z1)

REAL(EB), INTENT(OUT) :: Z0,Z1
REAL(EB) :: U1,U2,A

CALL RANDOM_NUMBER(U1)
CALL RANDOM_NUMBER(U2)
A = SQRT(-2._EB*LOG(U1))
Z0 = A*COS(TWOPI*U2)
Z1 = A*SIN(TWOPI*U2)

END SUBROUTINE BOX_MULLER


!> \brief Flux-limiting function

REAL(EB) FUNCTION F_B(B)
REAL(EB), INTENT(IN) :: B

IF (B<=0._EB) THEN
   F_B = 1._EB
ELSE
   F_B = (1._EB+B)**0.7_EB*LOG(1._EB+B)/B
ENDIF
END FUNCTION F_B


END MODULE MATH_FUNCTIONS



!> \brief Functions for physical quantities

MODULE PHYSICAL_FUNCTIONS

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
IMPLICIT NONE

CONTAINS


!> \brief Check if the species mass fractions are in bounds
!> \param ZZ_IN Array of mass fractions

LOGICAL FUNCTION IS_REALIZABLE(ZZ_IN)

REAL(EB), INTENT(IN) :: ZZ_IN(1:N_TRACKED_SPECIES)
REAL(EB), PARAMETER:: ZERO_MINUS=-EPSILON(0._FB),ONE_PLUS=1._EB+EPSILON(1._FB)

IF (ANY(ZZ_IN<ZERO_MINUS) .OR. SUM(ZZ_IN)>ONE_PLUS) THEN
   IS_REALIZABLE=.FALSE.
ELSE
   IS_REALIZABLE=.TRUE.
ENDIF

END FUNCTION IS_REALIZABLE


!> \brief Clip mass fractions between zero and one and redistribute clipped mass to most abundant species
!> \param ZZ_GET Array of mass fractions

SUBROUTINE GET_REALIZABLE_MF(ZZ_GET)

REAL(EB), INTENT(INOUT) :: ZZ_GET(1:N_TRACKED_SPECIES)
REAL(EB) :: SUM_OTHER_SPECIES
INTEGER :: N_ZZ_MAX

! clip mass fractions
ZZ_GET=MAX(0._EB,MIN(1._EB,ZZ_GET))

! absorb all error in most abundant species
N_ZZ_MAX = MAXLOC(ZZ_GET,1)

SUM_OTHER_SPECIES = SUM(ZZ_GET) - ZZ_GET(N_ZZ_MAX)
ZZ_GET(N_ZZ_MAX) = 1._EB - SUM_OTHER_SPECIES

END SUBROUTINE GET_REALIZABLE_MF


!> \brief Determine the mass fraction of a primitive species
!> \param Z_IN Array of lumped mass fractions
!> \param INDEX Index of desired primitive species
!> \param Y_OUT Mass fraction of desired primitive species

SUBROUTINE GET_MASS_FRACTION(Z_IN,INDEX,Y_OUT)

INTEGER, INTENT(IN) :: INDEX
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT

Y_OUT = DOT_PRODUCT(Z2Y(INDEX,1:N_TRACKED_SPECIES),Z_IN)
Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))

END SUBROUTINE GET_MASS_FRACTION


!> \brief Determine the mass fractions of all primitive species
!> \param Z_IN Array of lumped species mass fractions
!> \param Y_OUT Array of primitive species mass fractions

SUBROUTINE GET_MASS_FRACTION_ALL(Z_IN,Y_OUT)

REAL(EB), INTENT(IN)  :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT(1:N_SPECIES)
INTEGER :: I

DO I=1,N_SPECIES
   Y_OUT(I) = DOT_PRODUCT(Z2Y(I,1:N_TRACKED_SPECIES),Z_IN)
ENDDO

Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))

END SUBROUTINE GET_MASS_FRACTION_ALL


!> \brief Compute ratio of molecular weight of all gas components except the one specified to the mol. wgt. of the one specified
!> \param INDEX_IN Index of the gas species
!> \param MW_RATIO W_all/W_in
!> \param Y_IN Optional mass fraction of primitive species
!> \param Z_IN Optional mass fraction of lumped species

SUBROUTINE GET_MW_RATIO(INDEX_IN,MW_RATIO,Y_IN,Z_IN)

REAL(EB), INTENT(IN), OPTIONAL :: Y_IN(1:N_SPECIES), Z_IN(1:N_TRACKED_SPECIES)
INTEGER, INTENT(IN):: INDEX_IN
REAL(EB), INTENT(OUT) :: MW_RATIO
INTEGER:: NS

MW_RATIO = 0._EB
IF (PRESENT(Y_IN)) THEN
   IF (ABS(Y_IN(INDEX_IN)-1._EB) > TWO_EPSILON_EB) THEN
      DO NS=1,N_SPECIES
         IF (NS==INDEX_IN) CYCLE
         MW_RATIO = MW_RATIO + Y_IN(NS)/SPECIES(NS)%MW
      ENDDO
      IF (MW_RATIO<=TWO_EPSILON_EB) THEN
         MW_RATIO=SPECIES_MIXTURE(1)%MW
      ELSE
         MW_RATIO = (1._EB-Y_IN(INDEX_IN))/MW_RATIO
      ENDIF
   ELSE
      MW_RATIO=SPECIES_MIXTURE(1)%MW
   ENDIF
   MW_RATIO = MW_RATIO/SPECIES(INDEX_IN)%MW
ELSE
   IF (ABS(Z_IN(INDEX_IN)-1._EB) > TWO_EPSILON_EB) THEN
      DO NS=1,N_TRACKED_SPECIES
         IF (NS==INDEX_IN) CYCLE
         MW_RATIO = MW_RATIO + Z_IN(NS)/SPECIES_MIXTURE(NS)%MW
      ENDDO
      IF (MW_RATIO<=TWO_EPSILON_EB) THEN
         MW_RATIO=SPECIES_MIXTURE(1)%MW
      ELSE
         MW_RATIO = (1._EB-Z_IN(INDEX_IN))/MW_RATIO
      ENDIF
   ELSE
      MW_RATIO=SPECIES_MIXTURE(1)%MW
   ENDIF
   MW_RATIO = MW_RATIO/SPECIES_MIXTURE(INDEX_IN)%MW
ENDIF

END SUBROUTINE GET_MW_RATIO


SUBROUTINE GET_EQUIL_DATA(MW,TMP_L,PRES_IN,H_V,H_V_EFF,T_BOIL_EFF,X_EQ,H_V_LOWER,H_V_DATA,RAMP_INDEX,CONSTANT_H)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM
REAL(EB), INTENT(IN):: MW,TMP_L,PRES_IN
REAL(EB), INTENT(IN), OPTIONAL :: H_V_DATA(:),CONSTANT_H
REAL(EB), INTENT(INOUT):: T_BOIL_EFF
REAL(EB), INTENT(OUT):: H_V,H_V_EFF,X_EQ
INTEGER, INTENT(IN), OPTIONAL:: H_V_LOWER, RAMP_INDEX
REAL(EB):: DHOR

IF (PRESENT(RAMP_INDEX)) THEN
   H_V=EVALUATE_RAMP(TMP_L,0._EB,RAMP_INDEX)
   H_V_EFF=EVALUATE_RAMP(T_BOIL_EFF,0._EB,RAMP_INDEX)
   DHOR = H_V_EFF*MW/R0
   T_BOIL_EFF = MAX(0._EB,DHOR*T_BOIL_EFF/(DHOR-T_BOIL_EFF*LOG(PRES_IN/P_STP)+TWO_EPSILON_EB))
   H_V_EFF=EVALUATE_RAMP(T_BOIL_EFF,0._EB,RAMP_INDEX)
ELSEIF (PRESENT(CONSTANT_H)) THEN
   H_V=CONSTANT_H
   H_V_EFF=CONSTANT_H
   DHOR = H_V_EFF*MW/R0
   T_BOIL_EFF = MAX(0._EB,DHOR*T_BOIL_EFF/(DHOR-T_BOIL_EFF*LOG(PRES_IN/P_STP)+TWO_EPSILON_EB))
ELSE
   CALL INTERPOLATE1D_UNIFORM(H_V_LOWER,H_V_DATA,TMP_L,H_V)
   CALL INTERPOLATE1D_UNIFORM(H_V_LOWER,H_V_DATA,T_BOIL_EFF,H_V_EFF)
   DHOR = H_V_EFF*MW/R0
   T_BOIL_EFF = MAX(0._EB,DHOR*T_BOIL_EFF/(DHOR-T_BOIL_EFF*LOG(PRES_IN/P_STP)+TWO_EPSILON_EB))
   CALL INTERPOLATE1D_UNIFORM(H_V_LOWER,H_V_DATA,T_BOIL_EFF,H_V_EFF)
ENDIF
H_V_EFF = 0.5_EB*(H_V+H_V_EFF)
X_EQ = MIN(1._EB,EXP(H_V_EFF*MW/R0*(1._EB/T_BOIL_EFF-1._EB/TMP_L)))

END SUBROUTINE GET_EQUIL_DATA


!> \brief Compute volume fraction of a primitive species
!> \param Y_INDEX Index of primitive species
!> \param ZZ_ARRAY Array of lumped species
!> \param R_MIX R/W of species mixture

REAL(EB) FUNCTION GET_VOLUME_FRACTION(Y_INDEX,ZZ_ARRAY,R_MIX)

INTEGER, INTENT(IN) :: Y_INDEX
REAL(EB), INTENT(IN) :: R_MIX,ZZ_ARRAY(:)
REAL(EB) :: MASS_FRACTION,RCON

CALL GET_MASS_FRACTION(ZZ_ARRAY,Y_INDEX,MASS_FRACTION)
RCON = SPECIES(Y_INDEX)%RCON
GET_VOLUME_FRACTION = RCON*MASS_FRACTION/R_MIX

END FUNCTION GET_VOLUME_FRACTION


!> \brief Determine molecular weight of the gas mixture
!> \param Z_IN Array of lumped species mass fractions
!> \param MW_OUT Average molecular weight (g/mol)

SUBROUTINE GET_MOLECULAR_WEIGHT(Z_IN,MW_OUT)

REAL(EB), INTENT(IN)  :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: MW_OUT

MW_OUT =  1._EB/DOT_PRODUCT(MWR_Z,Z_IN)

END SUBROUTINE GET_MOLECULAR_WEIGHT


!> \brief Compute R/W for a gas mixture
!> \param Z_IN Array of lumped species mass fractions
!> \param RSUM_OUT R/W_avg

SUBROUTINE GET_SPECIFIC_GAS_CONSTANT(Z_IN,RSUM_OUT)

REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: RSUM_OUT

RSUM_OUT =  R0 * DOT_PRODUCT(MWR_Z,Z_IN)

END SUBROUTINE GET_SPECIFIC_GAS_CONSTANT


!> \brief Get specific heat of the gas mixture at a specified temperature
!> \param Z_IN Array of lumped species mass fractions
!> \param CP_OUT Specific heat of the mixture (J/kg/K)
!> \param TMPG Gas mixture temperature (K)

SUBROUTINE GET_SPECIFIC_HEAT(Z_IN,CP_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: CP_OUT

ITMP = MIN(5000,NINT(TMPG))
CP_OUT  = DOT_PRODUCT(CP_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)

END SUBROUTINE GET_SPECIFIC_HEAT


!> \brief Get sensible enthalpy of the gas mixture at a specified temperature
!> \param Z_IN Array of lumped species mass fractions
!> \param H_S_OUT Specific heat of the mixture (J/kg)
!> \param TMPG Gas mixture temperature (K)

SUBROUTINE GET_SENSIBLE_ENTHALPY(Z_IN,H_S_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: H_S_OUT

ITMP = MIN(4999,INT(TMPG))

H_S_OUT = DOT_PRODUCT(H_SENS_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)+(TMPG-REAL(ITMP,EB))* &
          DOT_PRODUCT(H_SENS_Z(ITMP+1,1:N_TRACKED_SPECIES)-H_SENS_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)

END SUBROUTINE GET_SENSIBLE_ENTHALPY


SUBROUTINE GET_GIBBS_FREE_ENERGY(G_F_OUT,NU_IN,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG,NU_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: G_F_OUT

ITMP = MIN(5000,NINT(TMPG))
G_F_OUT = 1.E6_EB*DOT_PRODUCT(G_F_Z(ITMP,1:N_TRACKED_SPECIES),NU_IN) ! convert from kJ/mol to J/kmol

END SUBROUTINE GET_GIBBS_FREE_ENERGY


!> \brief Get average specific heat of the gas mixture up to a specified temperature
!> \param Z_IN Array of lumped species mass fractions
!> \param CPBAR_OUT Average specific heat of the mixture (J/kg/K)
!> \param TMPG Gas mixture temperature (K)

SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT(Z_IN,CPBAR_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: CPBAR_OUT

ITMP = MIN(5000,NINT(TMPG))
CPBAR_OUT = DOT_PRODUCT(CPBAR_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)

END SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT


!> \brief Get sensible enthalpy of lumped species N at a specified temperature
!> \param N Index of lumped species
!> \param TMPG Gas mixture temperature (K)
!> \param H_S Sensible enthalpy of lumped species N (J/kg)

SUBROUTINE GET_SENSIBLE_ENTHALPY_Z(N,TMPG,H_S)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
INTEGER, INTENT(IN) :: N
REAL(EB), INTENT(OUT) :: H_S

ITMP = MIN(4999,INT(TMPG))
H_S = H_SENS_Z(ITMP,N)+(TMPG-REAL(ITMP,EB))*(H_SENS_Z(ITMP+1,N)-H_SENS_Z(ITMP,N))

END SUBROUTINE GET_SENSIBLE_ENTHALPY_Z


!> \brief Get average specific heat of lumped species N up to a specified temperature
!> \param N Index of lumped species
!> \param TMPG Gas mixture temperature (K)
!> \param CPBAR_OUT Average specific heat of lumped species N (J/kg/K)

SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT_Z(N,TMPG,CPBAR_OUT)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
INTEGER, INTENT(IN) :: N
REAL(EB), INTENT(OUT) :: CPBAR_OUT

ITMP = MIN(5000,NINT(TMPG))
CPBAR_OUT = CPBAR_Z(ITMP,N)

END SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT_Z


!> \brief Get thermal conductivity of gas mixture at a specified temperature
!> \param Z_IN Array of lumped species mass fractions
!> \param K_OUT Conductivity of gas mixture (W/m/K)
!> \param TMPG Gas mixture temperature (K)

SUBROUTINE GET_CONDUCTIVITY(Z_IN,K_OUT,TMPG)

REAL(EB), INTENT(IN) :: Z_IN(1:N_TRACKED_SPECIES), TMPG
REAL(EB), INTENT(OUT) :: K_OUT
INTEGER :: ITMP

ITMP = MIN(5000,NINT(TMPG))
K_OUT = DOT_PRODUCT(K_RSQMW_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)/DOT_PRODUCT(Z_IN,RSQ_MW_Z)

END SUBROUTINE GET_CONDUCTIVITY


SUBROUTINE GET_SOLID_CONDUCTIVITY(K_OUT,TMP_S,OPT_MATL_INDEX,OPT_SURF_INDEX,OPT_RHO_IN,OPT_I_IN,OPT_J_IN,OPT_K_IN)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER, INTENT(IN), OPTIONAL :: OPT_MATL_INDEX,OPT_SURF_INDEX,OPT_I_IN,OPT_J_IN,OPT_K_IN
REAL(EB), INTENT(IN), OPTIONAL :: OPT_RHO_IN(1:N_MATL)
REAL(EB), INTENT(IN) :: TMP_S
REAL(EB), INTENT(OUT) :: K_OUT
INTEGER :: N,NR,I_LOC,J_LOC,K_LOC
REAL(EB) :: VOLSUM
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

!! K_OUT = 0.2_EB; RETURN ! PMMA debug

K_OUT = 0._EB
IF (PRESENT(OPT_MATL_INDEX)) THEN
   ML => MATERIAL(OPT_MATL_INDEX)
   IF (ML%K_S>0._EB) THEN
      K_OUT = ML%K_S
   ELSE
      NR = -NINT(ML%K_S)
      K_OUT = EVALUATE_RAMP(TMP_S,0._EB,NR)
   ENDIF
ELSEIF (PRESENT(OPT_SURF_INDEX)) THEN
   ! See FDS Tech Guide, Eqs. (7.22)-(7.25)
   SF => SURFACE(OPT_SURF_INDEX)
   K_OUT = 0._EB
   VOLSUM = 0._EB
   DO N=1,SF%N_MATL
      IF (OPT_RHO_IN(N)<TWO_EPSILON_EB) CYCLE
      ML => MATERIAL(SF%MATL_INDEX(N))
      VOLSUM = VOLSUM + OPT_RHO_IN(N)/ML%RHO_S
      IF (ML%K_S>0._EB) THEN
         K_OUT = K_OUT + OPT_RHO_IN(N)*ML%K_S/ML%RHO_S
      ELSE
         NR = -NINT(ML%K_S)
         K_OUT = K_OUT + OPT_RHO_IN(N)*EVALUATE_RAMP(TMP_S,0._EB,NR)/ML%RHO_S
      ENDIF
   ENDDO
   IF (VOLSUM>TWO_EPSILON_EB) THEN
      K_OUT = K_OUT/VOLSUM
   ! ! debug
   ! ELSE
   !    print *,SF%ID,SF%N_MATL
   ENDIF
ENDIF
IF (K_OUT<TWO_EPSILON_EB) THEN
   IF (PRESENT(OPT_I_IN)) I_LOC = OPT_I_IN
   IF (PRESENT(OPT_J_IN)) J_LOC = OPT_J_IN
   IF (PRESENT(OPT_K_IN)) K_LOC = OPT_K_IN
   ! ! debug
   ! WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY K_OUT: ',K_OUT
   ! WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY VOLSUM: ',VOLSUM
   ! IF (PRESENT(OPT_I_IN)) WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY I: ',I_LOC
   ! IF (PRESENT(OPT_J_IN)) WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY J: ',J_LOC
   ! IF (PRESENT(OPT_K_IN)) WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY K: ',K_LOC
   K_OUT = 10000._EB
ENDIF

END SUBROUTINE GET_SOLID_CONDUCTIVITY


SUBROUTINE GET_SOLID_RHOCBAR(RHOCBAR_OUT,TMP_S,OPT_MATL_INDEX,OPT_SURF_INDEX,OPT_RHO_IN)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER, INTENT(IN), OPTIONAL :: OPT_MATL_INDEX,OPT_SURF_INDEX
REAL(EB), INTENT(IN), OPTIONAL :: OPT_RHO_IN(1:N_MATL)
REAL(EB), INTENT(IN) :: TMP_S
REAL(EB), INTENT(OUT) :: RHOCBAR_OUT
INTEGER :: N,NR
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

!! RHOCBAR_OUT = 2420000._EB; RETURN ! PMMA debug

RHOCBAR_OUT = 0._EB
IF (PRESENT(OPT_MATL_INDEX)) THEN
   ML => MATERIAL(OPT_MATL_INDEX)
   IF (ML%C_S>0._EB) THEN
      RHOCBAR_OUT = ML%RHO_S*ML%C_S
   ELSE
      NR = -NINT(ML%C_S)
      RHOCBAR_OUT = ML%RHO_S*EVALUATE_RAMP(TMP_S,0._EB,NR)
   ENDIF
ELSEIF (PRESENT(OPT_SURF_INDEX)) THEN
   ! See FDS Tech Guide, Eq. (7.24)
   SF => SURFACE(OPT_SURF_INDEX)
   DO N=1,SF%N_MATL
      IF (OPT_RHO_IN(N)<TWO_EPSILON_EB) CYCLE
      ML => MATERIAL(SF%MATL_INDEX(N))
      IF (ML%C_S>0._EB) THEN
         RHOCBAR_OUT = RHOCBAR_OUT + OPT_RHO_IN(N)*ML%C_S
      ELSE
         NR = -NINT(ML%C_S)
         RHOCBAR_OUT = RHOCBAR_OUT + OPT_RHO_IN(N)*EVALUATE_RAMP(TMP_S,0._EB,NR)
      ENDIF
   ENDDO
ENDIF
IF (RHOCBAR_OUT<=TWO_EPSILON_EB) RHOCBAR_OUT = 0.001_EB

END SUBROUTINE GET_SOLID_RHOCBAR


SUBROUTINE GET_SOLID_ABSORPTION_COEFFICIENT(KAPPA_OUT,SURF_INDEX,RHO_IN)

INTEGER, INTENT(IN) :: SURF_INDEX
REAL(EB), INTENT(IN) :: RHO_IN(1:N_MATL)
REAL(EB), INTENT(OUT) :: KAPPA_OUT
INTEGER :: N
REAL(EB) :: VOLSUM
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

KAPPA_OUT = 0._EB
VOLSUM = 0._EB
SF => SURFACE(SURF_INDEX)
DO N=1,SF%N_MATL
   IF (RHO_IN(N)<=TWO_EPSILON_EB) CYCLE
   ML => MATERIAL(SF%MATL_INDEX(N))
   VOLSUM = VOLSUM + RHO_IN(N)/ML%RHO_S
   KAPPA_OUT = KAPPA_OUT + RHO_IN(N)*ML%KAPPA_S/ML%RHO_S
ENDDO
IF (VOLSUM>0._EB) KAPPA_OUT = KAPPA_OUT/VOLSUM

END SUBROUTINE GET_SOLID_ABSORPTION_COEFFICIENT


SUBROUTINE GET_SOLID_REFRACTIVE_INDEX(INDEX_OUT,SURF_INDEX,RHO_IN)

INTEGER, INTENT(IN) :: SURF_INDEX
REAL(EB), INTENT(IN) :: RHO_IN(1:N_MATL)
REAL(EB), INTENT(OUT) :: INDEX_OUT
INTEGER :: N
REAL(EB) :: VOLSUM
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

INDEX_OUT = 0._EB
VOLSUM = 0._EB
SF => SURFACE(SURF_INDEX)
DO N=1,SF%N_MATL
   IF (RHO_IN(N)<=TWO_EPSILON_EB) CYCLE
   ML => MATERIAL(SF%MATL_INDEX(N))
   VOLSUM = VOLSUM + RHO_IN(N)/ML%RHO_S
   INDEX_OUT = INDEX_OUT + RHO_IN(N)*ML%REFRACTIVE_INDEX/ML%RHO_S
ENDDO
IF (VOLSUM>0._EB) INDEX_OUT = INDEX_OUT/VOLSUM

END SUBROUTINE GET_SOLID_REFRACTIVE_INDEX


!> \brief Get viscosity of gas mixture at a specified temperature
!> \param Z_IN Array of lumped species mass fractions
!> \param MU_OUT Viscosity of gas mixture (kg/m/s)
!> \param TMPG Gas mixture temperature (K)

SUBROUTINE GET_VISCOSITY(Z_IN,MU_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN)  :: TMPG,Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: MU_OUT

ITMP = MIN(5000,NINT(TMPG))
MU_OUT = DOT_PRODUCT(MU_RSQMW_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)/DOT_PRODUCT(Z_IN,RSQ_MW_Z)

END SUBROUTINE GET_VISCOSITY


SUBROUTINE GET_ENTHALPY(Z_IN,H_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES),DTMP
REAL(EB), INTENT(OUT) :: H_OUT

IF (TMPG>=5000._EB) THEN
   H_OUT = DOT_PRODUCT(CPBAR_Z(5000,1:N_TRACKED_SPECIES),Z_IN)*TMPG
ELSE
   ITMP = INT(TMPG)
   DTMP = TMPG-REAL(ITMP)
   H_OUT = DOT_PRODUCT(CPBAR_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)
   H_OUT = H_OUT+DTMP*(DOT_PRODUCT(CPBAR_Z(ITMP+1,1:N_TRACKED_SPECIES),Z_IN)-H_OUT)
   H_OUT = H_OUT*TMPG
ENDIF

END SUBROUTINE GET_ENTHALPY


REAL(EB) FUNCTION DRAG(RE,DRAG_LAW)

! drag coefficient

INTEGER, INTENT(IN) :: DRAG_LAW
REAL(EB), INTENT(IN) :: RE

SELECT CASE(DRAG_LAW)

   ! see J.P. Holman 7th Ed. Fig. 6-10
   CASE(SPHERE_DRAG)
      IF (RE<=TWO_EPSILON_EB) THEN
         DRAG = 100._EB
      ELSEIF (RE<=1._EB) THEN
         DRAG = 24._EB/RE
      ELSEIF (RE<1000._EB) THEN
         !!DRAG = 24._EB*(1._EB+0.15_EB*RE**0.687_EB)/RE ! see Crowe, Sommerfeld, Tsuji, 1998, Eq. (4.51)
         DRAG = 24._EB*(0.85_EB+0.15_EB*RE**0.687_EB)/RE ! matches Stokes drag at RE=1 (RJM)
      ELSEIF (RE>=1000._EB) THEN
         DRAG = 0.44_EB
      ENDIF

   ! see J.P. Holman 7th Ed. Fig. 6-9
   CASE(CYLINDER_DRAG)
      IF (RE<=1._EB) THEN
         DRAG = 10._EB/(RE**0.8_EB)
      ELSEIF (RE>1._EB .AND. RE<1000._EB) THEN
         DRAG = 10._EB*(0.6_EB+0.4_EB*RE**0.8_EB)/RE
      ELSEIF (RE>=1000._EB) THEN
         DRAG = 1._EB
      ENDIF

   CASE(USER_DRAG)
      DRAG = 1._EB ! PC%DRAG_COEFFICIENT set elsewhere

   CASE DEFAULT
      DRAG = 0._EB

END SELECT

END FUNCTION DRAG


!> \brief Compute the surface (mass or energy) density of a wall cell
!> \param NM Mesh number
!> \param MODE Indicator of returned quantity (0) kg/m2, (1) kg/m3, (2) kJ/m2, (3) kJ/m3
!> \param WALL_INDEX Index of the WALL cell
!> \param LAGRANGIAN_PARTICLE_INDEX Index of the particle
!> \param MATL_INDEX Index of the material component

REAL(EB) FUNCTION SURFACE_DENSITY(NM,MODE,WALL_INDEX,LAGRANGIAN_PARTICLE_INDEX,MATL_INDEX)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER, INTENT(IN) :: NM,MODE
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,LAGRANGIAN_PARTICLE_INDEX,MATL_INDEX
INTEGER :: I_GRAD=0,NWP,II2,N,NR
REAL(EB) :: WGT,R_S(0:NWP_MAX),EPUM
TYPE(MESH_TYPE), POINTER :: M=>NULL()
TYPE(ONE_D_M_AND_E_XFER_TYPE), POINTER :: ONE_D=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (MATERIAL_TYPE),POINTER :: ML=>NULL()

M => MESHES(NM)

IF (PRESENT(WALL_INDEX)) THEN
   ONE_D => M%WALL(WALL_INDEX)%ONE_D
   SF => SURFACE(M%WALL(WALL_INDEX)%SURF_INDEX)
ELSEIF (PRESENT(LAGRANGIAN_PARTICLE_INDEX)) THEN
   ONE_D => M%LAGRANGIAN_PARTICLE(LAGRANGIAN_PARTICLE_INDEX)%ONE_D
   SF => SURFACE(LAGRANGIAN_PARTICLE_CLASS(M%LAGRANGIAN_PARTICLE(LAGRANGIAN_PARTICLE_INDEX)%CLASS_INDEX)%SURF_INDEX)
ELSE
   SURFACE_DENSITY = 0
   RETURN
ENDIF

THERMALLY_THICK_IF: IF (SF%THERMAL_BC_INDEX/=THERMALLY_THICK) THEN

   SURFACE_DENSITY = 0._EB

ELSE THERMALLY_THICK_IF

   SELECT CASE(SF%GEOMETRY)
      CASE(SURF_CARTESIAN)    ; I_GRAD = 1
      CASE(SURF_CYLINDRICAL)  ; I_GRAD = 2
      CASE(SURF_SPHERICAL)    ; I_GRAD = 3
   END SELECT

   NWP = SUM(ONE_D%N_LAYER_CELLS)
   R_S(0:NWP) = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(0:NWP)

   SURFACE_DENSITY = 0._EB
   NUMBER_WALL_POINTS_LOOP: DO II2=1,NWP
      AREA_VOLUME_SELECT: SELECT CASE(MODE)
         CASE(0,2); WGT = (R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(REAL(I_GRAD,EB)*(SF%INNER_RADIUS+SF%THICKNESS)**(I_GRAD-1))
         CASE(1,3); WGT = (R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(SF%INNER_RADIUS+SF%THICKNESS)**I_GRAD
      END SELECT AREA_VOLUME_SELECT

      EPUM = 1._EB ! energy per unit mass
      IF (PRESENT(MATL_INDEX)) THEN
         ENERGY_SELECT_1: SELECT CASE(MODE)
            CASE(2,3)
               ML => MATERIAL(MATL_INDEX)
               IF (ML%C_S>0._EB) THEN
                  EPUM = ML%C_S * ONE_D%TMP(II2)
               ELSE
                  NR = -NINT(ML%C_S)
                  EPUM = EVALUATE_RAMP(ONE_D%TMP(II2),0._EB,NR) * ONE_D%TMP(II2)
               ENDIF
         END SELECT ENERGY_SELECT_1
         SURFACE_DENSITY = SURFACE_DENSITY + ONE_D%MATL_COMP(MATL_INDEX)%RHO(II2)*WGT*EPUM
      ELSE
         DO N=1,SF%N_MATL
            ENERGY_SELECT_2: SELECT CASE(MODE)
               CASE(2,3)
                  ML => MATERIAL(N)
                  IF (ML%C_S>0._EB) THEN
                     EPUM = ML%C_S * ONE_D%TMP(II2)
                  ELSE
                     NR = -NINT(ML%C_S)
                     EPUM = EVALUATE_RAMP(ONE_D%TMP(II2),0._EB,NR) * ONE_D%TMP(II2)
                  ENDIF
            END SELECT ENERGY_SELECT_2
            SURFACE_DENSITY = SURFACE_DENSITY + ONE_D%MATL_COMP(N)%RHO(II2)*WGT*EPUM
         ENDDO
      ENDIF

   ENDDO NUMBER_WALL_POINTS_LOOP

ENDIF THERMALLY_THICK_IF

END FUNCTION SURFACE_DENSITY


!> \brief Compute particle Cumulative Number Fraction (CNF) and Cumulative Volume Fraction (CVF)
!> \param DM Median particle diameter (m)
!> \param RR Array of particle radii (m)
!> \param CNF Cumulative Number Fraction
!> \param CVF Cumulative Volume Fraction
!> \param NPT Number of points in the distribution
!> \param GAMMA Parameter in the distribution function
!> \param SIGMA Parameter in the distribution function
!> \param DISTRIBUTION Type of distribution

SUBROUTINE PARTICLE_SIZE_DISTRIBUTION(DM,RR,CNF,CVF,NPT,GAMMA,SIGMA,DISTRIBUTION)

USE MATH_FUNCTIONS, ONLY: IERFC
CHARACTER(LABEL_LENGTH), INTENT(IN) :: DISTRIBUTION
REAL(EB), INTENT(IN) :: DM,GAMMA,SIGMA
INTEGER, INTENT(IN) :: NPT
REAL(EB) :: SUM1,SUM2,DD1,DI,ETRM,GFAC,SFAC,DMIN,X1
INTEGER  :: J
REAL(EB), INTENT(OUT) :: RR(0:NPT),CNF(0:NPT),CVF(0:NPT)

RR(0)  = 0._EB
CNF(0) = 0._EB
SUM1   = 0._EB
SUM2   = 0._EB

X1     = IERFC(2._EB*CNF_CUTOFF)
DMIN   = MAX(DM*EXP(-X1*SQRT(2._EB)*SIGMA),0._EB)
DD1    = (-LOG(CNF_CUTOFF)/LOG(2._EB))**(1._EB/GAMMA)*DM
DD1    = (DD1-DMIN)/REAL(NPT,EB)
GFAC   = LOG(2._EB)*GAMMA*DD1/(DM**GAMMA)
SFAC   = DD1/(SQRT(TWOPI)*SIGMA)

INTLOOP: DO J=1,NPT
   DI = DMIN + (J-0.5_EB)*DD1
   RR(J) = 0.5_EB*DI
   IF ((DI<=DM .OR. DISTRIBUTION=='LOGNORMAL') .AND. DISTRIBUTION/='ROSIN-RAMMLER') THEN
      ETRM = EXP(-(LOG(DI/DM))**2/(2._EB*SIGMA**2))
      SUM1 = SUM1 + (SFAC/DI**4)*ETRM
      SUM2 = SUM2 + (SFAC/DI)*ETRM
   ELSE
      ETRM = EXP(-LOG(2._EB)*(DI/DM)**GAMMA)
      SUM1 = SUM1 + GFAC*DI**(GAMMA-4._EB)*ETRM
      SUM2 = 1._EB - ETRM
   ENDIF
   CNF(J) = SUM1
   CVF(J) = SUM2
ENDDO INTLOOP

CNF = CNF/SUM1
CVF = CVF/SUM2

END SUBROUTINE PARTICLE_SIZE_DISTRIBUTION


SUBROUTINE SPRAY_ANGLE_DISTRIBUTION(LON,LAT,LON_CDF,LAT_CDF,BETA,MU,SPRAY_ANGLE,DISTRIBUTION_TYPE,NPT)

INTEGER,INTENT(IN) :: NPT
REAL(EB),INTENT(OUT) :: LON_CDF(0:NPT),LON(0:NPT),LAT(0:NPT),LAT_CDF(0:NPT,0:NPT)
REAL(EB),INTENT(IN) :: BETA,MU,SPRAY_ANGLE(2,2)
CHARACTER(LABEL_LENGTH),INTENT(IN) :: DISTRIBUTION_TYPE
INTEGER :: I,J
REAL(EB) :: dLON,dLAT,PDF(0:NPT,0:NPT),THETA_MAX,THETA_MIN

THETA_MAX=MAXVAL(SPRAY_ANGLE)
THETA_MIN=MINVAL(SPRAY_ANGLE)

dLAT=(THETA_MAX-THETA_MIN)/REAL(NPT)
dLON=2._EB*PI/REAL(NPT)
!Discretize latitude and longtitude
LAT=(/ (THETA_MIN+I*dLAT,I=0,NPT) /)
LON=(/ (0._EB+I*dLON, I=0,NPT) /)

! SPRAY_ANGLE May be different in X and Y directions
! i.e spray angle is dependent on longtitude
! SPRAY_AGLE_MIN and MAX form ellipses with semi-axes defined by
! SPRAY_ANGLE(1,1:2) and SPRAY_ANGLE(2,1:2) respectively

DO I=0,NPT
  THETA_MIN=SPRAY_ANGLE(1,1)*SPRAY_ANGLE(1,2)
  IF(THETA_MIN>0._EB) THEN
     THETA_MIN=THETA_MIN/SQRT((SPRAY_ANGLE(1,2)*COS(LON(I)))**2+(SPRAY_ANGLE(1,1)*SIN(LON(I)))**2)
  ENDIF
  THETA_MAX=SPRAY_ANGLE(2,1)*SPRAY_ANGLE(2,2)
  THETA_MAX=THETA_MAX/SQRT((SPRAY_ANGLE(2,2)*COS(LON(I)))**2+(SPRAY_ANGLE(2,1)*SIN(LON(I)))**2)
  SELECT CASE(DISTRIBUTION_TYPE)
   CASE("TRIANGLE")
      DO J=0,NPT
         IF(LAT(J)<THETA_MIN .OR. LAT(J)>THETA_MAX) THEN
           PDF(J,I)=0._EB
         ELSE
           IF(LON(I)<MU) THEN
              PDF(J,I)=2*(LAT(J)-THETA_MIN)/((THETA_MAX-THETA_MIN)*(THETA_MAX-MU))
           ELSE
              PDF(J,I)=2*(THETA_MAX-LAT(J))/((THETA_MAX-THETA_MIN)*(THETA_MAX-MU))
           ENDIF
         ENDIF
      ENDDO
   CASE("GAUSSIAN")
      DO J=0,NPT
        IF(LAT(J)<THETA_MIN .OR. LAT(J)>THETA_MAX) THEN
           PDF(J,I)=0._EB
        ELSE
           PDF(J,I)=exp(-BETA*((LAT(J)-MU)/(THETA_MAX-THETA_MIN))**2)
        ENDIF
      ENDDO
   CASE DEFAULT ! "UNIFORM"
       DO J=0,NPT
        IF(LAT(J)<THETA_MIN .OR. LAT(J)>THETA_MAX) THEN
           PDF(J,I)=0._EB
        ELSE
           PDF(J,I)=1._EB
        ENDIF
      ENDDO
   END SELECT
ENDDO


!
DO I=0,NPT
PDF(I,:)=PDF(I,:)*SIN(LAT(I))
ENDDO

LAT_CDF=0._EB
! Latitude distribution conditional on Longtitude
DO I=1,NPT
   LAT_CDF(I,:)=LAT_CDF(I-1,:)+0.5*(PDF(I,:)+PDF(I-1,:))*dLAT
ENDDO

! Marginal longtitude distribution
LON_CDF=0._EB
DO I=1,NPT
   LON_CDF(I)=LON_CDF(I-1)+0.5*(LAT_CDF(NPT,I-1)+LAT_CDF(NPT,I))*dLON
ENDDO

! Normalize marginal longtitude distribution
LON_CDF=LON_CDF/LON_CDF(NPT)
!Normalize conditional latitude distributions
DO I=1,NPT
   LAT_CDF(:,I)=LAT_CDF(:,I)/LAT_CDF(NPT,I)
ENDDO

END SUBROUTINE SPRAY_ANGLE_DISTRIBUTION


REAL(EB) FUNCTION FED(Y_IN,RSUM,FED_ACTIVITY)

! Returns the integrand of FED (Fractional Effective Dose) calculation.

REAL(EB), INTENT(IN) :: Y_IN(1:N_TRACKED_SPECIES),RSUM
INTEGER, INTENT(IN) :: FED_ACTIVITY
INTEGER  :: N
REAL(EB) :: Y_MF_INT, TMP_1
REAL(EB), DIMENSION(3) :: CO_FED_FAC
!                at rest           light work(default) heavy work
DATA CO_FED_FAC /0.70486250E-5_EB, 2.7641667E-5_EB,    8.2925E-5_EB/

! All equations from D.A. Purser, Sec. 2, Chap. 6, SFPE Handbook, 4th Ed.
! Note: Purser uses minutes, here dt is in seconds. Conversion at the end of the function.
! Total FED dose:
! FED_dose = (FED_LCO + FED_LCN + FED_LNOx + FLD_irr)*FED_VCO2 + FED_LO2;

FED = 0._EB

! Carbon monoxide (CO)
!          at rest    light work heavy work
! RMV_FED /8.5_EB,    25.0_EB,   50.0_EB /
! D_FED   /40.0_EB,   30.0_EB,   20.0_EB /
! RMV/D   /0.2125_EB, 0.8333_EB, 2.5_EB/
!
! FED_LCO = (3.317E-5 * (C_CO)^1.036 * RMV * (dt/60)) / D;
!   with RMV=25 [l/min], D=30 [%] COHb concentration at incapacitation and C_CO in ppm
!
IF (CO_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,CO_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(CO_INDEX)%RCON*Y_MF_INT*1.E6_EB/RSUM
   ! FED   = 2.764E-5_EB*TMP_1**(1.036_EB)
   FED   = CO_FED_FAC(FED_ACTIVITY)*TMP_1**(1.036_EB)
ENDIF

! Nitrogen oxides (NOx, here NO + NO2)
! FED_LNOx = C_NOx/1500 * (dt/60);
!   with C_NOx = C_NO + C_NO2, all in ppm
TMP_1 = 0._EB
IF (NO_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,NO_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(NO_INDEX)%RCON*Y_MF_INT/RSUM
ENDIF
IF (NO2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,NO2_INDEX,Y_MF_INT)
   TMP_1 = TMP_1 + SPECIES(NO2_INDEX)%RCON*Y_MF_INT/RSUM
ENDIF
IF (TMP_1 > 0._EB) FED = FED + TMP_1/0.001500_EB

! Cyanide
! FED_LCN = (exp(C_CN/43)/220 - 0.0045) * (dt/60);
!   with C_CN = C_HCN - C_NOx, all in ppm
IF (HCN_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,HCN_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(HCN_INDEX)%RCON*Y_MF_INT/RSUM - TMP_1
   IF (TMP_1 > 0._EB) FED = FED + (Exp(TMP_1/0.000043_EB)/220.0_EB-0.00454545_EB)
ENDIF

! Irritants
! FLD_irr = (C_HCl/F_HCl + C_HBr/F_HBr + C_HF/F_HF + C_SO2/F_SO2 + C_NO2/F_NO2 + C_C3H4O/F_C3H4O + C_CH2O/F_CH2O) * (dt/60);
!   all in ppm
TMP_1 = 0._EB
DO N=1,N_SPECIES
   IF (SPECIES(N)%FLD_LETHAL_DOSE > 0._EB) THEN
      Call GET_MASS_FRACTION(Y_IN,N,Y_MF_INT)
      TMP_1 = TMP_1 + SPECIES(N)%RCON*Y_MF_INT/RSUM / SPECIES(N)%FLD_LETHAL_DOSE
   ENDIF
ENDDO
FED = FED + TMP_1

! Carbon dioxide (CO2) induced hyperventilation:
! FED_VCO2 = exp(0.1903*C_CO2/1E4 + 2.0004)/7.1;
!   C_CO2 in ppm
IF (CO2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,CO2_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(CO2_INDEX)%RCON*Y_MF_INT/RSUM
   If ( TMP_1 > 0.0_EB ) FED = FED * Exp( 0.1903_EB*TMP_1*100.0_EB + 2.0004_EB )/7.1_EB
ENDIF

! Low oxygen (O2)
! FED_LO2 = 1/exp(8.13 - 0.54*(0.209 - C_O2/1E6)) * (dt/60);
!   C_O2 in ppm
IF (O2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,O2_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(O2_INDEX)%RCON*Y_MF_INT/RSUM
   IF ( TMP_1 < 0.20_EB ) FED = FED + 1.0_EB  / Exp(8.13_EB-0.54_EB*(20.9_EB-100.0_EB*TMP_1))
ENDIF

! Convert the FED integrand for minutes.
FED = FED / 60._EB

END FUNCTION FED


REAL(EB) FUNCTION FIC(Y_IN,RSUM)
! Returns FIC (Fractional Incapacitating Concentration)

REAL(EB), INTENT(IN) :: Y_IN(1:N_TRACKED_SPECIES),RSUM
REAL(EB) :: Y_MF_INT
INTEGER  :: N

FIC = 0._EB
DO N=1,N_SPECIES
   IF (SPECIES(N)%FIC_CONCENTRATION > 0._EB) THEN
      Call GET_MASS_FRACTION(Y_IN,N,Y_MF_INT)
      FIC = FIC + SPECIES(N)%RCON*Y_MF_INT/RSUM / SPECIES(N)%FIC_CONCENTRATION
   ENDIF
ENDDO

END FUNCTION FIC


REAL(EB) FUNCTION WATER_VAPOR_MASS_FRACTION(HUMIDITY,TEMP,PZONE)

! Compute the water vapor mass fraction given the relative humidity and temperature

REAL(EB), INTENT(IN) :: TEMP,HUMIDITY,PZONE
REAL(EB) :: X_SAT,DHOR_T,DHOR,P_RATIO,T_BOIL_EFF
REAL(EB),PARAMETER :: T_BOIL=373.15_EB,DHOR_T_B=4916.346083_EB

DHOR_T = (H_V_H2O(INT(TEMP))+(TEMP-REAL(INT(TEMP,EB)))*(H_V_H2O(INT(TEMP)+1)-H_V_H2O(INT(TEMP))))*MW_H2O/R0
P_RATIO = PZONE/P_STP
T_BOIL_EFF = MAX(0._EB,DHOR_T_B*T_BOIL/(DHOR_T_B-T_BOIL*LOG(P_RATIO)+TWO_EPSILON_EB))
DHOR = 0.5_EB*((H_V_H2O(INT(T_BOIL_EFF))+(T_BOIL_EFF-REAL(INT(T_BOIL_EFF,EB)))*&
       (H_V_H2O(INT(T_BOIL_EFF)+1)-H_V_H2O(INT(T_BOIL_EFF))))*MW_H2O/R0+DHOR_T)
X_SAT  = MIN(1._EB,EXP(DHOR*(1._EB/T_BOIL_EFF-1._EB/TEMP)))
WATER_VAPOR_MASS_FRACTION = HUMIDITY*0.01_EB*X_SAT/(MW_AIR/MW_H2O+(1._EB-MW_AIR/MW_H2O)*HUMIDITY*0.01_EB*X_SAT)

END FUNCTION WATER_VAPOR_MASS_FRACTION


REAL(EB) FUNCTION RELATIVE_HUMIDITY(Y_H2O,TEMP,PZONE)

! Compute the relative humidity given the water vapor mass fraction and temperature

REAL (EB), INTENT(IN) :: TEMP,Y_H2O,PZONE
REAL (EB) :: X_SAT,X_H2O,DHOR,DHOR_T,T_BOIL_EFF,P_RATIO
REAL(EB),PARAMETER :: T_BOIL=373.15_EB,DHOR_T_B=4916.346083_EB

P_RATIO = PZONE/P_STP
T_BOIL_EFF = MAX(0._EB,DHOR_T_B*T_BOIL/(DHOR_T_B-T_BOIL*LOG(P_RATIO)+TWO_EPSILON_EB))
IF (TEMP >= T_BOIL_EFF) THEN
   X_SAT = 1._EB
ELSE
   DHOR_T = (H_V_H2O(INT(TEMP))+(TEMP-REAL(INT(TEMP,EB)))*(H_V_H2O(INT(TEMP)+1)-H_V_H2O(INT(TEMP))))*MW_H2O/R0
   DHOR = 0.5_EB*((H_V_H2O(INT(T_BOIL_EFF))+(T_BOIL_EFF-REAL(INT(T_BOIL_EFF,EB)))*&
          (H_V_H2O(INT(T_BOIL_EFF)+1)-H_V_H2O(INT(T_BOIL_EFF))))*MW_H2O/R0+DHOR_T)
   X_SAT  = MIN(1._EB,EXP(DHOR*(1._EB/T_BOIL_EFF-1._EB/TEMP)))
ENDIF
X_H2O = Y_H2O*MW_AIR/(MW_H2O-Y_H2O*(MW_H2O-MW_AIR))
RELATIVE_HUMIDITY = 100._EB * X_H2O / X_SAT

END FUNCTION RELATIVE_HUMIDITY


REAL(EB) FUNCTION LES_FILTER_WIDTH_FUNCTION(DX,DY,DZ)
   REAL(EB), INTENT(IN):: DX,DY,DZ
   IF (TWO_D) THEN
      LES_FILTER_WIDTH_FUNCTION = SQRT(DX*DZ)
   ELSE
      LES_FILTER_WIDTH_FUNCTION = (DX*DY*DZ)**ONTH
   ENDIF
END FUNCTION LES_FILTER_WIDTH_FUNCTION


REAL(EB) FUNCTION GET_POTENTIAL_TEMPERATURE(TMP_IN,Z_IN)

USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP
REAL(EB), INTENT(IN) :: TMP_IN,Z_IN
REAL(EB) :: PP,DUMMY

IF (STRATIFICATION) THEN
   PP = EVALUATE_RAMP(Z_IN,DUMMY,I_RAMP_P0_Z)
ELSE
   PP = P_INF
ENDIF
GET_POTENTIAL_TEMPERATURE = TMP_IN*(1.E5_EB/PP)**GM1OG  ! GM1OG = (GAMMA-1)/GAMMA = R/CP

END FUNCTION GET_POTENTIAL_TEMPERATURE


SUBROUTINE MONIN_OBUKHOV_SIMILARITY(Z,Z_0,L,U_STAR,THETA_STAR,THETA_0,U,TMP)

REAL(EB), INTENT(IN) :: Z,Z_0,L,U_STAR,THETA_STAR,THETA_0
REAL(EB), INTENT(OUT) :: U,TMP
REAL(EB), PARAMETER :: P_REF=100000._EB,RHO_REF=1.2_EB
REAL(EB) :: PSI_M,PSI_H,THETA,KAPPA

KAPPA = VON_KARMAN_CONSTANT
CALL MONIN_OBUKHOV_STABILITY_CORRECTIONS(PSI_M,PSI_H,Z,L)
U = (U_STAR/KAPPA)*(LOG(Z/Z_0)-PSI_M)
THETA = THETA_0 + (THETA_STAR/KAPPA)*(LOG(Z/Z_0)-PSI_H)
TMP = THETA*(P_REF/(P_REF-RHO_REF*GRAV*Z))**(-0.285_EB)

END SUBROUTINE MONIN_OBUKHOV_SIMILARITY


SUBROUTINE MONIN_OBUKHOV_STABILITY_CORRECTIONS(PSI_M,PSI_H,Z,L)

! Reference: A.J. Dyer. A review of flux profile relationships. Boundary-Layer Meteorology, 7:363–372, 1974.

REAL(EB), INTENT(IN) :: Z,L
REAL(EB), INTENT(OUT) :: PSI_M,PSI_H
REAL(EB) :: ZETA

IF (L>=0._EB) THEN
   ! stable boundary layer
   PSI_M = -5._EB*Z/L
   PSI_H = PSI_M
ELSE
   ! unstable boundary layer
   ZETA = (1._EB-16._EB*Z/L)**0.25_EB
   PSI_M = 2._EB*LOG(0.5_EB*(1._EB+ZETA)) + LOG(0.5_EB*(1._EB+ZETA**2)) - 2._EB*ATAN(ZETA) + 0.5_EB*PI
   PSI_H = 2._EB*LOG(0.5_EB*(1._EB+ZETA**2))
ENDIF

END SUBROUTINE MONIN_OBUKHOV_STABILITY_CORRECTIONS


!> \brief Computes the mass and heat transfer coeffiicents for a liquid droplet in air based on the selected EVAP_MODEL on MISC
!> \param H_MASS The mass transfer coefficient (m2/s)
!> \param H_HEAT The dropelt heat transfer coefficient (W/m2/K)
!> \param D_AIR Diffusivity in air of the droplet species in the gas cell with the droplet (m2/s)
!> \param K_AIR Conductivity in the gas cell with the droplet (W/m/k)
!> \param CP_AIR Specific heat in the gas cell with the droplet (J/kg/K)
!> \param RHO_AIR Density in the gas cell with the droplet (kg/m3)
!> \param LENGTH Length scale (m)
!> \param Y_DROP Equilibrium vapor fraction for the current droplet temperature
!> \param Y_GAS Mass fraction of vapor in the gas
!> \param B_NUMBER B number for the droplet
!> \param NU_FAC_GAS Constant factor used in computing  the NUsselt number
!> \param SH_FAC_GAS Constant factor used in computing  the Sherwood number
!> \param RE_L Renyolds number
!> \param TMP_FILM Film temperature for the droplet (K)
!> \param ZZ_GET Tracked species mass fractions in the gas cell with the droplet
!> \param Z_INDEX Droplet species index in ZZ

SUBROUTINE DROPLET_H_MASS_H_HEAT_GAS(H_MASS,H_HEAT,D_AIR,K_AIR,CP_AIR,RHO_AIR,LENGTH,Y_DROP,Y_GAS,B_NUMBER,NU_FAC_GAS,SH_FAC_GAS, &
                                     RE_L,TMP_FILM,ZZ_GET,Z_INDEX)
USE MATH_FUNCTIONS, ONLY: F_B
REAL(EB), INTENT(IN) :: D_AIR,CP_AIR,K_AIR,RHO_AIR,LENGTH,Y_DROP,Y_GAS,NU_FAC_GAS,SH_FAC_GAS,RE_L,TMP_FILM, &
                        ZZ_GET(1:N_TRACKED_SPECIES)
INTEGER, INTENT(IN) :: Z_INDEX
REAL(EB), INTENT(INOUT) :: B_NUMBER
REAL(EB), INTENT(OUT) :: H_MASS,H_HEAT
REAL(EB) :: NUSSELT,SHERWOOD,LEWIS,THETA,C_GAS_DROP,C_GAS_AIR,ZZ_GET2(1:N_TRACKED_SPECIES)

SELECT CASE (EVAP_MODEL)
   CASE(-1) ! Ranz Marshall
      NUSSELT  = 2._EB + NU_FAC_GAS*SQRT(RE_L)
      H_HEAT   = NUSSELT*K_AIR/LENGTH
      IF (Y_DROP <= Y_GAS) THEN
         H_MASS   = 0._EB
      ELSE
         SHERWOOD = 2._EB + SH_FAC_GAS*SQRT(RE_L)
         H_MASS   = SHERWOOD*D_AIR/LENGTH
      ENDIF
   CASE(0) ! Sazhin M0, Eq 106 + 109 with B_T=B_M. This is the default model.
      IF (Y_DROP <= Y_GAS) THEN
         NUSSELT  = 2._EB + NU_FAC_GAS*SQRT(RE_L)
         H_HEAT   = NUSSELT*K_AIR/LENGTH
         H_MASS   = 0._EB
      ELSE
         NUSSELT  = ( 2._EB + NU_FAC_GAS*SQRT(RE_L) )*LOG(1._EB+B_NUMBER)/B_NUMBER
         H_HEAT   = NUSSELT*K_AIR/LENGTH
         SHERWOOD = ( 2._EB + SH_FAC_GAS*SQRT(RE_L) )*LOG(1._EB+B_NUMBER)/(Y_DROP-Y_GAS)
         H_MASS   = SHERWOOD*D_AIR/LENGTH
         ! above we save a divide and multiply of B_NUMBER
         ! the full model corresponding to Sazhin (108) and (109) would be
         ! SH = SH_0 * LOG(1+B_M)/B_M
         ! H_MASS = SH * D/L * B_M/(Y_D-Y_G)
      ENDIF
   CASE(1) ! Sazhin M1, Eq 106 + 109 with eq 102.
      IF (Y_DROP <= Y_GAS) THEN
         NUSSELT  = 2._EB + NU_FAC_GAS*SQRT(RE_L)
         H_HEAT   = NUSSELT*K_AIR/LENGTH
         H_MASS   = 0._EB
      ELSE
         SHERWOOD = ( 2._EB + SH_FAC_GAS*SQRT(RE_L) )*LOG(1._EB+B_NUMBER)/(Y_DROP-Y_GAS)
         H_MASS   = SHERWOOD*D_AIR/LENGTH
         LEWIS    = K_AIR / (RHO_AIR * D_AIR * CP_AIR)
         ZZ_GET2(1:N_TRACKED_SPECIES) = 0._EB
         ZZ_GET2(Z_INDEX) = 1._EB
         CALL GET_SPECIFIC_HEAT(ZZ_GET2,C_GAS_DROP,TMP_FILM)
         CALL GET_SPECIFIC_HEAT(ZZ_GET,C_GAS_AIR,TMP_FILM)
         THETA = C_GAS_DROP/C_GAS_AIR/LEWIS
         B_NUMBER = (1._EB+B_NUMBER)**THETA-1._EB
         NUSSELT  = ( 2._EB + NU_FAC_GAS*SQRT(RE_L) )*LOG(1._EB+B_NUMBER)/B_NUMBER
         H_HEAT   = NUSSELT*K_AIR/LENGTH
      ENDIF
   CASE(2) ! Sazhin M2, Eq 116 and 117 with eq 106, 109, and 102.
      IF (Y_DROP <= Y_GAS) THEN
         NUSSELT  = 2._EB + NU_FAC_GAS*SQRT(RE_L)
         H_HEAT   = NUSSELT*K_AIR/LENGTH
         H_MASS   = 0._EB
      ELSE
         SHERWOOD = ( 2._EB + SH_FAC_GAS*SQRT(RE_L) )*LOG(1._EB+B_NUMBER)/((Y_DROP-Y_GAS)*F_B(B_NUMBER))
         H_MASS   = SHERWOOD*D_AIR/LENGTH
         LEWIS    = K_AIR / (RHO_AIR * D_AIR * CP_AIR)
         ZZ_GET2(1:N_TRACKED_SPECIES) = 0._EB
         ZZ_GET2(Z_INDEX) = 1._EB
         CALL GET_SPECIFIC_HEAT(ZZ_GET2,C_GAS_DROP,TMP_FILM)
         CALL GET_SPECIFIC_HEAT(ZZ_GET,C_GAS_AIR,TMP_FILM)
         THETA = C_GAS_DROP/C_GAS_AIR/LEWIS
         B_NUMBER = (1._EB+B_NUMBER)**THETA-1._EB
         NUSSELT  = ( 2._EB + NU_FAC_GAS*SQRT(RE_L) )*LOG(1._EB+B_NUMBER)/(B_NUMBER*F_B(B_NUMBER))
         H_HEAT   = NUSSELT*K_AIR/LENGTH
      ENDIF
END SELECT

END SUBROUTINE DROPLET_H_MASS_H_HEAT_GAS

END MODULE PHYSICAL_FUNCTIONS



!> \brief Coordinate transformation functions

MODULE TRAN

USE PRECISION_PARAMETERS
IMPLICIT NONE
TYPE TRAN_TYPE
   REAL(EB), POINTER, DIMENSION(:,:) :: C1=>NULL(),C2=>NULL(),C3=>NULL(),CCSTORE=>NULL(),PCSTORE=>NULL()
   INTEGER, POINTER, DIMENSION(:,:) :: IDERIVSTORE=>NULL()
   INTEGER NOC(3),ITRAN(3),NOCMAX
END TYPE TRAN_TYPE
TYPE (TRAN_TYPE), ALLOCATABLE, TARGET, DIMENSION(:) :: TRANS


CONTAINS


!> \brief Transformation function that maps uniform numerical coordinate to non-uniform physical coordinate
!> \param X Coordinate on uniform numerical grid (m)
!> \param IC Coordinate index, 1, 2, or 3
!> \param NM Mesh number

REAL(EB) FUNCTION G(X,IC,NM)

REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T=>NULL()

T => TRANS(NM)

N = T%NOC(IC)
IF (N==0) THEN
   G = X
   RETURN
ENDIF

SELECT CASE(T%ITRAN(IC))
   CASE(1)
      G = 0._EB
      DO I=1,N+1
         G = G + T%C1(I,IC)*X**I
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      G = T%C2(II-1,IC) + T%C3(II,IC)*(X-T%C1(II-1,IC))
END SELECT

END FUNCTION G


!> \brief First derivative of the transformation function that maps uniform numerical coordinate to non-uniform physical coordinate
!> \param X Coordinate on uniform numerical grid (m)
!> \param IC Coordinate index, 1, 2, or 3
!> \param NM Mesh number

REAL(EB) FUNCTION GP(X,IC,NM)

REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T=>NULL()

T => TRANS(NM)
N =  T%NOC(IC)
IF (N==0) THEN
   GP = 1._EB
   RETURN
ENDIF

SELECT CASE(T%ITRAN(IC))
   CASE(1)
      GP = 0._EB
      DO I=1,N+1
         GP = GP + I*T%C1(I,IC)*X**(I-1)
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GP = T%C3(II,IC)
   CASE DEFAULT
      GP = 0._EB
END SELECT

END FUNCTION GP


!> \brief Inverse of the transformation function that maps uniform numerical coordinate to non-uniform physical coordinate
!> \param Z Coordinate on non-uniform physical grid (m)
!> \param IC Coordinate index, 1, 2, or 3
!> \param NM Mesh number

REAL(EB) FUNCTION GINV(Z,IC,NM)

REAL(EB) :: GF
INTEGER :: N,IT,II,I
REAL(EB), INTENT(IN) :: Z
INTEGER, INTENT(IN)  :: IC,NM
TYPE (TRAN_TYPE), POINTER :: T=>NULL()

T => TRANS(NM)
GINV = Z
N = T%NOC(IC)
IF (N==0) RETURN

SELECT CASE(T%ITRAN(IC))
   CASE(1)
      LOOP1: DO IT=1,10
         GF = G(GINV,IC,NM)-Z
         IF (ABS(GF)<0.00001_EB) EXIT LOOP1
         GINV = GINV - GF/GP(GINV,IC,NM)
      ENDDO LOOP1
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (Z<=T%C2(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GINV = T%C1(II-1,IC) + (Z-T%C2(II-1,IC))/T%C3(II,IC)
END SELECT

END FUNCTION GINV


!> \brief Determine the continuous and discrete form of the cell indices of a given point
!> \param X x-coordinate of the point (m)
!> \param Y y-coordinate of the point (m)
!> \param Z z-coordinate of the point (m)
!> \param NM Mesh number
!> \param XI Real number indicating the cell index of the point in the x direction
!> \param YJ Real number indicating the cell index of the point in the y direction
!> \param ZK Real number indicating the cell index of the point in the z direction
!> \param I Cell index of the point in the x direction
!> \param J Cell index of the point in the y direction
!> \param K Cell index of the point in the z direction
!> \details For example, if the mesh is the unit cube with 10 cm cells, for the point (X,Y,Z)=(0.23,0.46,0.66),
!> GET_IJK would return (XI,JK,ZK)=(2.3,4.6,6.6) and (I,J,K)=(3,5,7).

SUBROUTINE GET_IJK(X,Y,Z,NM,XI,YJ,ZK,I,J,K)

USE MESH_VARIABLES
REAL(EB),INTENT(IN) :: X,Y,Z
INTEGER,INTENT(IN) :: NM
REAL(EB), INTENT(OUT) :: XI,YJ,ZK
INTEGER, INTENT(OUT) :: I,J,K

XI = MESHES(NM)%CELLSI(  MIN( MESHES(NM)%CELLSI_HI, MAX(MESHES(NM)%CELLSI_LO,FLOOR((X-MESHES(NM)%XS)*MESHES(NM)%RDXINT)) )  )
YJ = MESHES(NM)%CELLSJ(  MIN( MESHES(NM)%CELLSJ_HI, MAX(MESHES(NM)%CELLSJ_LO,FLOOR((Y-MESHES(NM)%YS)*MESHES(NM)%RDYINT)) )  )
ZK = MESHES(NM)%CELLSK(  MIN( MESHES(NM)%CELLSK_HI, MAX(MESHES(NM)%CELLSK_LO,FLOOR((Z-MESHES(NM)%ZS)*MESHES(NM)%RDZINT)) )  )
I = FLOOR(XI+1._EB)
J = FLOOR(YJ+1._EB)
K = FLOOR(ZK+1._EB)

END SUBROUTINE GET_IJK


END MODULE TRAN



!> \brief Module for various OpenMP functions

MODULE OPENMP

USE GLOBAL_CONSTANTS, ONLY : OPENMP_AVAILABLE_THREADS, OPENMP_USED_THREADS, OPENMP_USER_SET_THREADS, USE_OPENMP
!$ USE OMP_LIB
IMPLICIT NONE
PUBLIC OPENMP_INIT, OPENMP_SET_THREADS, OPENMP_PRINT_STATUS

CONTAINS

!> \brief Set the control flag USE_OPENMP if OpenMP is used.

SUBROUTINE OPENMP_INIT

!$OMP PARALLEL
!$OMP MASTER
!$ IF (OMP_GET_NUM_THREADS() /= 0) THEN
!$    USE_OPENMP = .TRUE.
!$    OPENMP_AVAILABLE_THREADS = OMP_GET_NUM_THREADS()
!$ ENDIF
!$OMP END MASTER
!$OMP BARRIER
!$OMP END PARALLEL

END SUBROUTINE OPENMP_INIT


!> \brief Change the number of OpenMP threads if set by the user in the input file.

SUBROUTINE OPENMP_SET_THREADS

!$IF (OPENMP_USER_SET_THREADS .EQV. .TRUE.) THEN
!$  IF (OPENMP_AVAILABLE_THREADS .NE. OPENMP_USED_THREADS) THEN
!$      CALL OMP_SET_NUM_THREADS(OPENMP_USED_THREADS)
!$  END IF
!$ELSE
!$  OPENMP_USED_THREADS = OPENMP_AVAILABLE_THREADS
!$  CALL OMP_SET_NUM_THREADS(OPENMP_USED_THREADS)
!$END IF

END SUBROUTINE OPENMP_SET_THREADS


!> \brief Write OpenMP status to standard error.

SUBROUTINE OPENMP_PRINT_STATUS
  USE GLOBAL_CONSTANTS, ONLY : LU_ERR, MYID, N_MPI_PROCESSES, VERBOSE
  INTEGER :: THREAD_ID

  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(THREAD_ID)

  THREAD_ID = 0
  !$ THREAD_ID = OMP_GET_THREAD_NUM()

  !$OMP CRITICAL
  IF (USE_OPENMP .AND. OPENMP_USED_THREADS>1 .AND. VERBOSE) WRITE(LU_ERR,91) " OpenMP thread ",THREAD_ID," of ",&
     OPENMP_USED_THREADS-1," assigned to MPI process ",MYID," of ",N_MPI_PROCESSES-1
  IF (.NOT.USE_OPENMP .AND. VERBOSE) WRITE(LU_ERR,92) " MPI process ",MYID," of ",N_MPI_PROCESSES-1
  !$OMP END CRITICAL

  !$OMP END PARALLEL

  91 FORMAT(A,I3,A,I3,A,I6,A,I6)
  92 FORMAT(A,I6,A,I6)

END SUBROUTINE OPENMP_PRINT_STATUS

END MODULE OPENMP



MODULE MISC_FUNCTIONS

USE PRECISION_PARAMETERS
IMPLICIT NONE

CONTAINS


SUBROUTINE WRITE_SUMMARY_INFO(LU)

USE GLOBAL_CONSTANTS
USE COMP_FUNCTIONS, ONLY : GET_DATE
USE MPI
CHARACTER(LABEL_LENGTH) :: DATE
INTEGER, INTENT(IN) :: LU
INTEGER :: MPIVERSION,MPISUBVERSION,MPILIBLENGTH,IERR
CHARACTER(LEN=MPI_MAX_LIBRARY_VERSION_STRING) :: MPILIBVERSION

CALL GET_DATE(DATE)

WRITE(LU,'(/A/)')      ' Fire Dynamics Simulator'
WRITE(LU,'(A,A)')      ' Current Date     : ',TRIM(DATE)
WRITE(LU,'(A,A)')      ' Revision         : ',TRIM(GITHASH_PP)
WRITE(LU,'(A,A)')      ' Revision Date    : ',TRIM(GITDATE_PP)
#ifdef COMPVER_PP
WRITE(LU,'(A,A)')      ' Compiler         : ',TRIM(COMPVER_PP)
#endif
WRITE(LU,'(A,A/)')     ' Compilation Date : ',TRIM(BUILDDATE_PP)
                      WRITE(LU,'(A,I6)')  ' MPI Enabled;    Number of MPI Processes:  ',N_MPI_PROCESSES
IF (.NOT. USE_OPENMP) WRITE(LU,'(A)')     ' OpenMP Disabled'
IF (USE_OPENMP)       WRITE(LU,'(A,I6)')  ' OpenMP Enabled; Number of OpenMP Threads: ',OPENMP_AVAILABLE_THREADS

CALL MPI_GET_LIBRARY_VERSION(MPILIBVERSION,MPILIBLENGTH,IERR)
CALL MPI_GET_VERSION(MPIVERSION,MPISUBVERSION,IERR)
WRITE(LU,'(/A,I1,A,I1)') ' MPI version: ',MPIVERSION,'.',MPISUBVERSION
WRITE(LU,'(A,A)') ' MPI library version: ',TRIM(MPILIBVERSION)

END SUBROUTINE WRITE_SUMMARY_INFO


!> \brief Finds the device or control function assoicated with an input
!> \param NAME Namelist ID
!> \param CTRL_ID String containing name of control function.
!> \param DEVC_ID String containing name of device.
!> \param DEVICE_INDEX Integer index locating DEVC_ID in the array of devices.
!> \param CONTROL_INDEX Integer index locating CTRL_ID in the array of control functions.
!> \param INPUT_INDEX The current count of inputs of type NAME.

SUBROUTINE SEARCH_CONTROLLER(NAME,CTRL_ID,DEVC_ID,DEVICE_INDEX,CONTROL_INDEX,INPUT_INDEX)

USE DEVICE_VARIABLES, ONLY: DEVICE,N_DEVC
USE CONTROL_VARIABLES, ONLY: CONTROL,N_CTRL
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
CHARACTER(*), INTENT (IN) :: NAME,CTRL_ID,DEVC_ID
INTEGER :: I, DEVICE_INDEX,CONTROL_INDEX
INTEGER , INTENT(IN) :: INPUT_INDEX

! There cannot be both a device and controller for any given entity

IF (DEVC_ID /= 'null' .AND. CTRL_ID /='null') THEN
   WRITE(MESSAGE,'(A,A,1X,I3,A)')  'ERROR: ',TRIM(NAME),INPUT_INDEX,' has both a device (DEVC) and a control (CTRL) specified'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF

! Search for device

IF (DEVC_ID /= 'null') THEN
   DO I=1,N_DEVC
      IF (DEVICE(I)%ID==DEVC_ID) THEN
         DEVICE_INDEX = I
         RETURN
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(A,A,A)')  'ERROR: DEVICE ',TRIM(DEVC_ID),' does not exist'
   CALL SHUTDOWN(MESSAGE)
ENDIF

! Search for controller

IF (CTRL_ID /= 'null') THEN
   DO I=1,N_CTRL
      IF (CONTROL(I)%ID==CTRL_ID) THEN
         CONTROL_INDEX = I
         RETURN
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(A,A,A)')  'ERROR: CONTROL ',TRIM(CTRL_ID),' does not exist'
   CALL SHUTDOWN(MESSAGE)
ENDIF

END SUBROUTINE SEARCH_CONTROLLER


!> \brief Returns true if current MPI process MYID controls mesh NM or its neighbors.
!> \param NM Mesh number.

LOGICAL FUNCTION PROCESS_MESH_NEIGHBORHOOD(NM)

USE MESH_VARIABLES
USE GLOBAL_CONSTANTS, ONLY: MYID,PROCESS,EVAC_PROCESS
INTEGER, INTENT(IN) :: NM
INTEGER :: N

PROCESS_MESH_NEIGHBORHOOD = .FALSE.
DO N=1,MESHES(NM)%N_NEIGHBORING_MESHES
   IF (MYID==PROCESS(MESHES(NM)%NEIGHBORING_MESH(N))) PROCESS_MESH_NEIGHBORHOOD = .TRUE.
ENDDO
IF (MYID==PROCESS(NM) .OR. MYID==EVAC_PROCESS) PROCESS_MESH_NEIGHBORHOOD = .TRUE.

END FUNCTION PROCESS_MESH_NEIGHBORHOOD


END MODULE MISC_FUNCTIONS
